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Scowstcr  intrusion,  also  synonymously  referred  to  in  some  literature  as  saltwater 
intrusion,  has  become  an  increasingly  dominant  constraint  on  the  utilization  of 
gmirndwoter  in  coastal  aquifers.  Climate  change  is  a natural  factor  that  hos  been 
identified  to  have  serious  impacts  on  the  hydrologic  regime  and  water  resources.  When 
coupled  with  exteosive  freshwater  withdrawal  through  pumping,  climate  change  could 
pose  a serious  threat  to  coastal  freshwater  resources. 

The  research  presented  b this  dissertation  has  three  basic  components:  generation 
of  climate  change  impact  scenarios,  simulation  of  sea  level  rise  sensitive  variable-density 
groundwater  flow,  and  opiimiznlion  of  freshwater  withdrawal  from  coastal  oquifers.  The 
climate  change  scenario  generation  module  was  developed  using  the  spatial  model bg 
capahiiiiies  of  raster  based  geographical  information  systems  (GIS).  The  module 
specifically  models  sea  level  rise  and  generates  new  sea  levels,  shoreline  translations,  and 
impoci  maps. 


T^e  vanable-densily  groundwater  flow  simulation  module  simulates  the  influence 
of  sea  level  rise  in  altering  the  groiuidiniter  flow  regime.  A scenario  of  climate  change 
was  selected  and  applied.  Extensive  sensitivity  analyses  were  performed  to  test  the 
significance  of  various  parameters. 

The  optimization  model  is  based  on  the  concept  of  genetic  algorithms.  The  use  of 
evolutionary  algorithms  to  solve  optimization  problems  in  water  resources  is  a fairly 
recent  approach.  The  capability  of  evolutionary  algorithms  to  handle  nomlinear  problems 
has  created  interest  in  the  application  of  these  algorithms  in  the  solution  of  simulation* 
optimization  problems. 

The  research  nrethodology  followed  a systematic  approach  thal  involved 
investigation  of  climate  change  impacts,  simulation  of  ffeshwaier*sallwatcr  flow,  and 
optimization  of  freewater  withdrawals  subject  to  the  various  constrainls  considered.  A 
threc'dimensional  variable  density  groundwater  flow  model  and  a cmas*sectiona[  model 
were  developed  to  simulate  seawater  Intrusion. 

The  primary  objective  of  the  research  was  to  provide  an  integrated  approach  to  the 
investigation  of  seawater  intrusion  in  coastal  aquifers  and  to  develop  a bcucr 
understanding  of  the  processes  involved.  The  results  obtained  ftom  the  test  cases 
considered  clearly  demonstrdte  the  significance  of  sea  level  rise  as  a factor  to  be 
considered  bt  the  development  of  management  strategies  for  coastal  aquifers. 


CHAPTER  I 
INTRODUCTION 

Water  is  0 preciom  natural  resource.  It  is  esiimalcd  that  99.4%  of  the 
available  water  resources  of  the  earth  (1.4x109  km’)  are  surface  water.  Groundwater 
accounts  only  for  0.6%  (9x10^  km’)  of  the  total.  From  the  total  amount  of  surface  water, 
most  of  it  (97%)  is  in  the  form  of  salt  water  in  oceans  and  inland  seas.  Fresh  surface 
water  accounts  for  only  2%oflhelolal  volume  of  water  (Bear  etal..  1999). 

Ifsalt  water  is  excluded,  the  proportion  of  surface  and  sub-surface  freshwater 
resources  becomes  78%  and  22%.  respectively.  Most  of  the  freshwater  is  in  the  form  of 
ice  locked  in  the  icecaps  of  the  polar  regions.  This  amount  is  estimated  to  be  77%  of  the 
total  volume  of  surface  freshwater.  The  fresh  surface-water  resources  that  are  readily 
usable  for  human  consumption  are  a very  small  proportion,  0J%  in  lakes  and  0.(XI3%  in 
streams,  whereas  the  groundwater  is  comparatively  a very  large  22%.  Thus,  because 
water  shortages  are  unavoidable  in  Ihc  foreseeable  Future  (Gleick.  1993),  it  would  be  very 
prudent  to  conclude  that  the  future  of  mankind  is  dependent  on  a sustainable  exploitation 
of  groundwater. 

In  order  to  manage  properly  Iresh  groundwater  resources  in  coastal  areas,  it  is 
neceasary  to  investigate  the  possibility  of  intrusion  of  salt  water  into  coastal  aquifers. 
Saltwater  imrusion  is  described  as  a mass  iransport  ofsaline  waters  into  rones  oflhc 
aquifer  that  were  previously  occupied  by  Iresh  waIcr(Stewar1,  1999).  This degradalioit 
of  Ihe  water  quality  is  eharacleriaed  by  increases  In  salinity  levels  that  exceed  the 


acceptable  valuer  for  potable  oragricultuniJ  uses.  Tbta  problem  requires  a greater 
aliemion  because  about  70%  of  the  world  population  lives  in  coastal  plains. 


Global  climate  change  is  one  of  the  moat  pressing  environmental  issues  we  are 
facing  today.  The  Intergovernmental  Panel  on  Climate  Change  (IPCC)  ha.s  undertaken  a 
comprehensive  assessment  of  the  scientilic  evidence  about  global  climate  change.  The 
results  ofihe  IPCC  osscssmenis  were  published  In  a series  of  reporu  since  1990. 
According  to  the  reports,  substantial  changes  in  the  environment  have  occurred  in  the 
past  cenluries  as  a result  of  human  activities  associated  with  odvonces  in  science  and 
technology  (IPCC,  1995a,  2001a).  There  is  also  an  overwhelming  body  of  research 
conducted  by  climate  researchers  from  around  the  world  suggesting  that  human  aclivities 
have  ccmribuled  slgnilicanily  to  climate  change  (Kondratyev  and  Cracknell.  I99S,  and 
Storch  and  Floser,  1999). 

Through  the  years,  various  agrieultunil  and  industrial  practices  have  led  to 
increased  atmospheric  concentrations  of  several  greenhouse  gases,  including  carbon 
dioxide,  methane,  nitrous  oxide,  and  chloronuorocarbons  (CPC)  in  the  lower  part  of  Ihc 
aimosphcrc.  The  basic  heat-trapping  property  of  these  greenhouse  gases  is  essentially 
undisputed,  although  there  is  considerable  scicniific  uncertainty  about  exactly  how  and 
when  the  earth's  climate  will  respond  u>  enhanced  greenhouse  gases  in  the  liiture. 

Climate  Change 

The  importance  of  studying  climate  chonge  with  the  aim  of  establishing  abetter 
underslanding  of  the  climatic  process  and  determining  the  comributions  ofnaiural  causes 
and  anthropogenic  effecls  has  been  emphasiaed  in  the  wake  of  the  growing  concern  about 
the  possible  effecis  of  CCb.  chloronuorocarbons  (CFCs),  and  olhcr  greenhouse  gases. 

The  physical  and  chemical  factors  ihal  offeci  climate  and  bring  about  change  have  been 


under  Invesiigaiion.  The  mejorilyofcllmalological  research  in  recent  years  has  been 
conducted  with  the  aim  of  predicting  the  impacts  of  human  activities  on  the  climate. 
However,  the  diflicultics  involved  in  studying  the  efTecls  orhuntan  activities  directly 
have  necessitated  the  modeling  of  climate  change  to  be  based  primarily  on  predicting  the 
clintale  change  due  to  natural  factors  and  then  introducing  anthropogenic  effects 
(Kondratyev  and  Cracknell,  1998). 

There  arc  two  main  modeling  approaches  currently  in  use  to  forecast  the  irnpact  of 
climate  change.  The  first  one  involves  working  out  analogues  from  past  climatic  and 
hydrometric  records.  The  second,  which  is  more  commonly  used,  employs  mathematical 
simulations  of  climate  known  as  general  circulation  models  or  global  climate  models 
(GCMs). 

Climatic  analogues  are  developed  by  chorssing  warmer  than  average  years  from  the 
historical  record  and  identifying  what  conditions  were  like  to  gel  an  idea  of  the  likely 
effects  of  global  warming.  The  limitation  of  this  method  Is  that  the  relatively  short  period 
of  climatic  data  does  not  represent  conditions  that  may  come  with  future  climare  change. 
For  example,  unprecedented  long-term  increases  in  temperature,  particularly  in 
wintertime,  may  not  be  in  past  records 

The  GCMs  simulate  atmospheric  circuladon,  the  energy  exchanges,  and  other 
important  land/oceaiValmosphere  interactions.  However,  they  cannot  as  yet  model  well 
other  small'Scalc  processes  such  as  biological  processes,  precipitation,  and  cloud  cover. 
These  processes  have  significant  effects  on  water  resources.  The  GCMs  in  current  use 
project  climate  ox'er  several  decades  to  more  than  a century  ("climate"  being  based  on 


avci^gc  weather  conditions  over  10  years  or  more)  but  give  only  large-scale  predictions, 
not  the  regional  ones  needed  for  planning. 

Results  Irom  the  models  used  by  climatologists  generally  agree  that  the  temperature 
increase,  as  a global  annual  average,  might  be  from  PC  to  4°C  by  the  year  2100.  They 
also  agree  that  the  effect  would  he  grealer  in  the  higher  latitudes,  especially  during  the 
winter  months  and  over  large  land  mosses.  The  warmer  temperature  would  trigger  oUicr 
changes,  such  os  a change  in  global  precipitation  patterns,  a decrease  in  snow  and  ice 
coverage,  and  a rise  in  sea  levels.  GCMs  simulabng  a clbnatc  that  is  based  on  a doubling 
of  CO:  suggest  a global  mean  increase  in  precipitation  and  evaporation  ofbetween  3 and 
15%.  Yet  the  more  uscliil  informalion.  i.e..  the  location  and  timing  of  these  changes,  is 
still  uncerlain. 

Climate  Change  DeiectioD  and  Modeling 

The  climate  system  is  characterized  by  a complex  system  of  interacting  internal  and 
external  components.  The  internal  inieracdve  components  include  the  atmosphere,  the 
oceans,  sea  ice,  the  land  and  its  features  (including  vegetation,  albedo,  biomass,  and 
ecosystems),  snow  cover,  land  Ice  (Including  Ihc  semi-permanent  ice  sheets  of  Antaielica 
and  Greenland  and  glaciers),  and  hydrology  (including  rivers,  lakes,  and  surface  and 
subsurface  water)  (IPCC,  1995a).  The  components  of  the  climate  system  that  ore 
customarily  regarded  os  external  to  the  system  include  the  Sun  and  its  output,  the  Earth's 
rotation,  Sun-Earth  geometry,  and  the  slowly  changing  orbit,  the  physical  components  of 
the  Earth  system  such  as  the  distribuiiott  of  land  and  ocean,  geographic  features  on  the 
land,  ocean  bottom  topography  and  basin  configurations,  and  the  mass  and  baste 
composition  of  the  atmosphere  and  ocean  (IPCC.  1996a).  Figure  1.1  is  a schematic 
represenlalicm  of  the  components  of  l)ie  global  climate  system. 


Insurance  and  other  Rnancial  services. 


Impacts  of  Climate  Change  on  Water  Resources 

Global  warming  due  lo  greenhouse  gases  and  anthropogenic  climate  change  is 
likely  10  have  slgnificam  impacts  on  the  hydrological  cycle  lIPCe,  1995a).  Water 
resoumes  of  the  world,  which  are  already  stressed  in  some  areas,  are  highly  vulnerable  to 
climaie  change  (IPCC,  2001b).  The  major  impacis  afciiniaie  change  on  waler  resources 
and  the  hydrological  cycle  include  the  rdlowing; 

• Increased  competition  for  waler  supply  between  various  consumptive  and  non- 
consumplive  uses, 

• Variability  m precipitation  and  evaporation, 

• Variability  in  aquifer  recharge. 

• Increase  in  sea  level  rise,  and 

• Loss  of  freshwater  resources  and  increase  in  saltwater  intrusion  in  coastal  aquifers. 
Prccipilation  variability.  Climate  models  in  current  use  are  still  unable  lo  make 

precise  regional  prediclions.  This  inability  of  the  models  and  the  inherenl  complexity  of 
Ihe  hydrological  cycle  have  made  better  prediction  a diflicult  task.  However,  ii  is 
understood  that  a change  in  precipiiaiion  may  alTcci  surface  wetness,  re  Declivity,  and 
vegetation,  which  then  alTeci  evapotranspimlion  and  cloud  formation,  which  in  nini 
a^eci  precipitattoQ.  Meanwhile,  the  hydrological  system  is  also  responding  to  other 
human  activities  such  as  deforestation,  urbanization,  and  ihe  over'Use  of  water  supplies 
loriliralion  variability.  Infiltration  is  one  of  Ihe  major  processes  of  the 
hydrological  cycle.  These  changes  in  seasonal  palletiis  may  afTecl  ihe  regional 
dislribution  of  both  ground  and  surface-waler  supplies.  Thmugh  inhltralion,  aquifers  gel 
rechurge,  and  the  vadose  zone  of  the  subsurface  retains  moisture  necessary  for 
agncullure.  especially  rain-fed  agriculture.  When  Ihe  precipitation  pattern  is  affceledby 
climate  change,  the  processes  of  infilinlion  and  recharge  also  will  be  affected.  Long 


spells  of  dry  seasons  result  in  Ihe  lowering  of  the  water  table  and  decrease  in  yield  of 
production  wells  used  for  water  supply  oragricultuie.  Several  models  also  suggest  that 
pmcipiuiion  will  become  more  intense  in  some  parts  ofthe  world,  and  os  a result,  there 
will  be  increases  in  the  magnitudes  of  surface  runolTand  floods. 

Increase  In  sea  level  rise.  Global  warming  has  been  attributed  to  the  thermal 
expansion  ofthe  ocean  water  and  the  melting  of  glacial  deposits  in  the  polar  regions. 
These  two  factors  are  the  major  conlributots  to  global  sea  level  rise,  which  is  expected  to 
affect  the  livelihood  of  people  living  in  coastal  areas.  Inhabitants  of  small  islands  might 
face  the  worst  consequences  of  all  as  a result  of  land  submergence  due  to  increased  sea 
level.  Incieasc  in  sea  level  also  affects  the  coastal  dynamics  of  low-lying  areas.  The 
processes  of  beach  erosion  and  deposition  will  be  altered  as  a result  of  sea  level  rise. 

Increase  la  seawater  lotnisioD.  The  increase  in  sea  level  and  the  subsequent 
submergence  of  coastal  areas  is  also  expected  to  result  in  the  loss  of  coastal  freshwater 
resources.  The  upstream  movement  of  saltwater  in  estuaries  will  affccl  freshwater  well 
Helds  upriver.  Coastal  aquifers,  which  are  already  under  sircss  in  many  areas,  maybe 
subject  to  Increased  levels  of  salinity  due  to  seawater  Intrusion. 

Sea  I.evel  Rise 

There  is  a general  consensus  about  the  issue  of  average  global  sea  level  rise.  There 
is  a substantial  wealth  of  research  findings  confirming  that  the  average  sea  level  has  been 
rising  globally  at  an  average  rate  of  about  to  2 mm  annually  for  at  least  the  last  century  or 
so  (Peltier  and  Tushingham.  1989;  Trupin  and  Wahr,  1990;  Douglas,  1991),  and  probably 
at  a much  smaller  rale  for  the  previous  several  millennia  (Flemming,  1978;  Flemming 
and  Webb.  I9gfi;  Kearney  and  Stevenson,  1991;  Shennan  and  Woodworth,  1992; 
Varekamp  el  al.,  1992).  Due  to  a number  of  comribuiing  factor?,  primarily  global 


warming,  moal  researchers  plausibly  argue  ihal  Ihe  rale  al  which  the  sea  level  will  nse 
globally  in  the  Inlure  will  be  much  greaier  than  the  present  rule. 

Sea  I.evel  Rise  Components 

Sea  level  rise  has  two  components.  Fiist,  what  is  known  as  eustatic  sea  level  rise  is 
the  result  of  global  warming  due  to  natural  and  anthropogenic  causes.  The  rise  in  global 
temperature  induces  thermal  expansion  of  Ihe  ocean  water,  melting  of  mountain  giaciers, 
and  possibly  the  melting  of  continenlal  ice  sheets  such  as  those  of  Greenland.  The 
second  component  is  Ihe  reluiive  sea  level  rise,  which  is  allribuled  to  the  vertical 
movement  of  the  land.  The  dynamic  nature  of  the  continents  has  resulted  in  continental 
movement  over  long  periods  of  time.  This  continental  movement  could  be  either  a rise  or 
iail  in  Ihe  land  surface  that  results  in  a relative  rise  or  fall  of  sea  level  with  respect  to  the 
land  surface.  Considering  a localized  scale,  relative  sea  level  rise  could  also  be  brought 
about  by  land  subsidence  as  a result  of  pumping  groundw  ater. 

RstimalloD  of  Sea  Level  Rise 

There  Is  a general  consensus  among  reaeaiehers  that  the  thermal  eapansion  lhat 
occurs  al  all  ocean  temperatures  is  the  major  contributing  factor  for  sea  level  rise  in  the 
20“  and  21”  centuries  (IPCC,  2001a).  The  IPCC  report  also  slates  that  the  change  in 
salinity  has  hod  a sign!  fleam  impact  on  Ihe  local  density  and  local  sea  level  but  little 
cITeci  on  Ihe  global  average  sea  level  change. 

The  analysis  of  observational  data  at  different  lime  scales  has  been  very 
instrumental  in  assessing  the  impact  of  climate  change  on  sea  level  rise.  On  Ihe 
geological  lime  scale,  evidence  from  coral  reefs,  oxygen  isotopes,  and  other  records  has 
demonsrrated  sufficiently  the  existence  of  a globally  coherent  connection  between  sea 
level  and  climate.  This  has  enabled  researchers  to  make  valid  speculations  about  an 


imminent  rise  in  global  sea  level  of  a few  meters  due  to  the  long  term  global  warming  of 
several  degrees  Celsius  (Tooley.  1993). 

A1  the  secular  lime  scale,’  tide-gauge  data  have  been  the  primary  sources  for  direct 
observationsof relative  sea  level  change  (Warrick,  1993).  On  the  basis  oflatc  Holocene 
sea  level  indicators  and  other  data,  GomiU  (1990)  concludes  thal  the  rate  of  global  sea 
level  rise  over  the  last  1 00  years  was  of  the  order  of  1 .0  lo  2.0  mm/yr.  Geological 
evidence  of  the  past  10,000  lo  20.000  years  indicates  that  major  temporal  and  spatial 
changes  occur  in  relotlve  sea  level  rise.  Observed  seo  level  ehonges  near  former  sites  of 
glaciation  are  brought  about  by  glacio-isostaiic  effect.^  On  the  other  hand,  changes 
observed  at  tectonically  stable  localities  far  (ram  the  former  ice  sheets  approximate  the 
global  averagesea  level  change  (IPCC.  2001a).  The  IPCC  (1996b)  estimates  dial  sea 
level  rise  will  be  In  the  range  of  20  lo  86  cm  from  1990  to  2100  with  a best  esllmaie  of49 
cm  (Warrick  et  al..  1996).  Analysis  of  data  from  sites  that  have  been  subject  lo  Ihc 
glacio-isosiatic  effect  indicate  that  the  ocean  volume  may  hove  increased  lo  add  2.5  to  3.5 
m to  the  global  average  sea  level  over  the  past  6,000  years  (IPCC,  2001a).  The 
geological  sea  level  esiimales  may  also  include  a component  (iom  thermal  expansion  and 
glacier  mass  changes.  The  contribution  of  ice  sheets  to  sen  level  rise  in  the  20th  and  21*’ 
centuries  in  response  to  earlier  climate  change  is  estimated  lo  be  limn  0.0  to  0.5  tnm/yr 
(IPCC.  2001a).  This  has  10  be  added  to  the  effect  ofihe  20’"  century  and  future  climate 
change,  AOGCM^  experiments  for  the  IS92a  scenario,  including  sulfate  aerosols. 


conduciul  for  Ihe  hundred-year  period  trom  1990  lo  2090  show  ihsi  ihe  globsl-svernge 
sea  level  rise  from  thermal  expansion  is  in  the  range  of  0.09  id  037  m.  The  results 
further  indicate  that  there  will  be  an  acceleration  through  the  21 century.  Expansion  for 
2040  to  2090  is  greater  than  for  1990  to  2040  by  a factor  of  1 .4  to  2. 1 (IPCC,  2001a). 

The  most  recent  assessment  for  Ihe  complete  range  of  AOGCMs  and  SRES^ 
scenarios,  including  uncertainties  in  land-ice  changes,  permafrost  changes,  and  sediment 
deposition,  projects  the  global  average  sea  level  rise  lobe  from  0.09  to  0.08  m over  the 
1990 10  2 100  estimate  with  a central  value  of  0.48  m (IPCC,  2001a).  The  central  value 
^ves  on  average  rale  of  2.2  lo4.4  times  the  rate  over  the  20th  century.  The  earlier 
assessment  conducted  using  time-dependent  sulfate  aerosols,  all  of  the  IS92  scenarios 
with  uncertainties,  and  a simple  model  with  climate  sensitivities  of  I.S  to  estimated 
the  global  average  sea  level  rise  to  be  in  the  range  of  0.13  to  0.94  m (Warrick  et  al., 

1996). 

Groundwater  flow  In  Coastal  Ar^uifers 

Coastal  aquifers  constirute  an  important  source  of  water  in  many  parts  of  the  world. 
In  addition  to  being  the  primary  sources  of  water,  the  proximity  of  coastal  aquifers  to  the 
seacoasi  puts  them  in  a unique  position.  The  geology  of  coastal  aquifers  and  their  unique 
relationship  to  the  marine  environment  make  the  investigation  of  groundwater  flow  in 
coastal  aquifers  a technically  challenging  task.  The  llow  of  groundwater  in  coastal 
aquifers  is  chaTacierized  by  variation  in  density,  which  becomes  more  pronounced  closer 
10  Ihe  coast.  The  existence  of  this  variation  in  density  is  brought  about  by  the  presence  of 
dissolved  solids  in  the  groundwater. 


[cd  10  the  IPCC  Special  itepon  on 


The  tolal  dissolved  solids  (TDS)  are  used  as  a measure  of  coucenlsaiion  and  density 
vanolion  in  groundwater.  The  most  important  of  all  the  dissolved  solids  in  terras  of 
chaiacieriaing  groundwater  quality  is  chloride  {Cf ).  In  order  to  characterize  the  type  of 
groundwater,  classiticalions  are  based  on  the  concentratiort  of  TDS  and  chloride.  Ttvo 


different  classification  tables  arc  provided  for  the  purpose  of  understanding.  The  first  one 
uses  the  TDS  as  a measure  of  concentration  while  the  second  one  uses  Cl  (see  Tables  l-l 
and  1*2).  Table  I-l  is  based  on  Fetter  (2001)  and  Table  1-2  is  based  on  Stuyzfand 
(1993). 


e 1-1.  Classification  ofwaterhased  on  total  dissolved  solids. 


Table  1-2.  Classification  of  fresh  and  saline  groundwater  on  the  basis  of  chloride 


Fresh-bracliish 

Braeltish 

Brackish-saline 

Hyperhaline  or  brine 


30-150 
150-300 
300-1,000 
1,000-  10,000 
10,000  - 20,000 
a2aooo 


In  addition  to  their  unique  geology  and  proximity  to  the  marine  environment, 
coastal  aquifers  are  also  susceptible  to  environmental  stresses,  anthropogenic  climate 
change,  and  human  activities.  Climate  change,  for  example,  alters  the  hydrologic 
characteristics  of  the  area  and  results  in  changes  in  precipitation,  evaporation,  and 
evapolranspiralion.  Sea  level  rise  and  changes  in  the  tide  gage  records  also  would  be 
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expected,  which  could  rexull  in  shoreline  retreat  and  flooding.  Excessive  pumping  of 
freshwater  aggravates  the  problem  hy  exacerbating  land  submergence  and  destabilizing 
the  balance  of  the  saltwater-freshwater  interface. 

There  arc  a number  of  factors  that  need  to  be  considered  on  different  scales  while 
trying  to  understand  the  flow  of  groundwater  in  coosutl  ac|uifere.  A thorough  analysis 
requires  investigating  the  hydrogeology  of  the  aquifer  using  available  geophysical 
methods,  analyzing  iheconsiUuentsoflhe  flow  In  lertnsofthe  dissolved  solids, 
identi  lying  the  various  boundary  conditions,  and  applying  a suitable  variable  density  flow 
model.  Figure  1 -2  (modified  from  Oude  Essink,  2000)  details  the  factors  that  need  to  be 
considered  in  the  investigation  of  groundwater  flow  in  coastal  aquifers. 

There  are  numerous  complexities  involved  in  modeling  density  dependent  flow  in 
porous  media.  The  complexities  range  from  the  inherent  properties  of  the  porous  medium 
to  the  mathematical  formulation  of  the  process  and  the  ensuing  difficulties  in  ohtatning 
analytical  solutions.  In  order  to  be  able  to  model  groundwater  flow,  a number  of  issues 
Involved  in  the  modeling  process  have  to  be  considered.  The  processes  involved  include 

• Undetslanding  the  hydrogeology  of  an  area, 

• Formulating  a well  posed  problem,  and 

• Selecting  eHicient  solution  techniques. 

tlydrogeology.  Proper  understanding  of  the  hydrogeology  of  an  area  to  be  modeled 
is  D crucial  factor  in  the  modeling  process.  Knowledge  of  the  inherent  hydraulic 
properties  of  the  porous  medium  such  as  porosity,  hydraulic  conductivity,  and 
transmissivity  and  geologic  factors  such  as  lithology.  straKgraphy,  and  the  geochemistry 
of  the  layer  greatly  facilitate  the  modeling  process  by  providing  the  necessary 
informaiian  about  the  subsurface  processes  that  govern  the  flow. 


13 


Figure  1-2.  Factor  affecting  ( 


Problem  rormuladon.  Forniulaiinge  well  posed  problem  greatly  enhances  the 

spatial  scales,  and  deciding  surTicienl  dimensions  to  appropriately  represent  the  flow  are 
some  of  the  issues  addressed  by  a well  posed  problem  The  processes  involved  in  the 
modeling  of  density  dependent  flow  such  as  advection,  difltision,  and  dispersion  and  the 
choice  of  appropriate  stale  and  boundary  conditions  constitute  the  major  pans  of  a well 
posed  problem. 

Solution  techniques.  Analytical  solutions  to  groundwater  flow  problems  are  very 
diRieulI  to  obtain  owing  lo  the  heterogeneity  of  the  medium  and  non-linearity  ofihe 
di^ereniial  equations.  This  problem  gels  even  more  complicated  when  the  flow  is  a 
density  dependent  flow.  As  a result,  for  most  two-dimensional  and  three-dimensional 
problems,  numerical  solutions  are  the  only  practical  way  lo  solve  the  problem.  The 
numerical  melhods  employed  can  uliUze  finite-element  or  finile-difrerence  discretization. 
Also,  the  iterative  processes  involved  in  the  solution  require  selection  of  an  efndeni 
numerical  algorithm.  These  numerical  algorithms  perform  the  iterative  temporal  and 
spatial  processes  on  the  basis  of  a specified  set  of  convcigenee  criteria.  The  spatial  and 
temporal  discretization  plays  a very  imponanl  role  in  the  accuracy  ofihe  result  and 
numerical  stability  ofihe  solution.  One  has  lo  keep  in  mind  dial  higher  spatial  and 
temporal  discretization  involves  higher  computing  and  slorage  capabilities,  and.  as  a 
resulq  faster  and  larger  memory  computers. 

Application  of  Geographic  Information  Systems  in  Water  Resources 
Over  the  past  couple  ofdecades.  the  held  of  Geographical  Information  Systems 
(GIS)  has  gained  an  increasing  applicabillly  in  ihe  held  ofhydrologic  modeling.  GIS  has 
proven  to  be  an  invaluable  tool  in  enhancing  the  spatial  dimensions  ofhydrologic  models. 


There  are  two  basic  dale  rormaia  in  CIS,  Raster  and  Vector.  Raster  CIS  works  on  the 
basis  of  cells  defined  by  rows  and  coiumns.  These  cells  contain  the  spatial  location  and 
refer  to  the  attribute  information  in  the  database  that  contains  the  layer  information. 
Vector  GIS,  on  the  other  hand,  uses  points,  iincs,  and  polygons  to  define  the  layer 
information. 

Raster  CIS  has  many  features  that  make  it  more  suitable  for  perrorming  hydrologic 
analyses.  These  advantages  are  as  follows: 

Speed.  Raster  GIS  analyses  lake  eousiderably  less  computer  processor  time  os 
compared  to  Vector  CIS  analyses. 

Adaptability.  Hydrologic  modeis  have  beneilted  from  remoteiy  sensed  data.  The 
dale  format  in  lemoieiy  sensed  data  such  as  sateilile  imagery  is  similar  to  the  data  format 
used  In  Raster  CIS,  Ihus  making  Raster  GIS  more  easily  adaptable  in  the  attempt  to 
couple  remote  sensing  and  CIS  wilh  hydrological  models. 

Suitability  to  represent  coatiDUOus  features.  The  spatial  disuibution  of 
hydrolo^c  parameters  is  not  defined  by  a clear-cut  boundary.  It  is  rather  continuous  in 
nature  and  crosses  boundaries  and  fealures.  This  unique  property  is  best  represented  by 
cell-based  data  formats  that  represent  the  spatial  gradual  change  in  data  attributes. 

Display  resolution.  The  resolution  of  computer  monitors  is  measured  in  terms  of 
pixels  that  are  also  cell  based.  Ihus  making  Raster  suitable  for  high  ^aiial  resolution 
analysis  and  display. 

Earlier  applications  of  GIS  ulIMzed  the  display  capabilily  of  Ihe  system.  In  other 
words,  model  resuils  were  Imported  Inlo  Ihe  database  and  displayed.  Ihus  enabling  one  to 
seethe  spatial  distribution  of  Ihcieaulls  of  the  modeling  process  across  Ihe  area  being 


16 

modeled.  Lalec  on.  derived  parameters  such  as  slope  and  aspect  were  first  computed  as 
pan  of  the  CIS  analysis,  and  then  the  results  were  used  in  hydrologic  models.  The  final 
analysis  was  again  displayed  using  the  advanced  cartographic  abilities  of  CIS. 

Advances  in  the  geo-processing  and  analytical  capabilities  of  the  CIS  software 
resulted  in  coupling  hydrological  models  with  CIS.  This  was  achieved  either  by  coupling 
the  hydrological  model  os  pan  of  a CIS  analysis  package  or  by  integrating  the  CIS 
module  within  the  framework  of  the  larger  hydrological  model.  This  requires  substantial 
programming  and  testing  In  order  to  ensure  a seamless  coupling  between  the  superior 
computational  capabilities  of  mathematical  models  and  the  spatial  modeling  capabtlilies 
of  CIS. 

The  advantages  of  this  method  rely  on  the  fact  that  the  CIS  component  can  be  used 
05  a pre-  and  post-processor  of  the  model.  As  a pre-processor,  layers  of  data  that  are 
spatially  referenced  to  the  same  geographical  location  are  created,  thus  minimizing  the 
need  for  an  external  data  set  lobe  called  during  execution  of  the  model.  As  a post- 
processor, the  GIS  component  imports  the  results  of  the  mathematical  model  analysis  and 
displays  it  and  also  performs  spatial  analysis  in  order  to  generate  derived  sets  of  dau. 
Spatially  Distributed  Modeling  in  GIS 

Most  hydrologic  models,  although  computationally  robust,  generally  experience 
difficuUies  in  addressing  the  issue  of  spatial  variability  of  the  processes  involved.  This 
problem  is  not  very  pronounced  in  lumped  models,  in  which  there  is  no  concern  about  the 
spatio]  variability  of  model  pararoeters.  On  the  other  hand,  distributed  models  rely  on  the 
accurate  representation  of  the  spaual  distribution  of  model  parameters  and  hydrologic 
response  units. 
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In  disiribulcd  modeling,  the  calchmeni  is  subdivided  into  smaller  units  of 
homogeneous  properties  and  the  contribution  of  each  smaller  unit,  referred  to  as  a 
Hydrologic  Response  Unit  (HRU),  is  eventually  accounted  for  to  simulate  the  response  of 
the  whole  caichmenl.  This  approach  comes  at  the  price  of  large  amounts  of  data  that 
need  to  be  collected,  stored,  processed,  and  presented.  IncToascd  spatial  variabiliry 
magnifies  the  scale  of  this  problem  significantly.  This  calls  for  a method  that  reduces  the 
size  while  maintaining  the  spatial  discretization  and  precision  required  for  a mote 
accurate  result. 

GIS  also  helps  integrate  a very  diverse  array  of  spatially  roferenced  data  through 
the  common  denominator,  which  is  their  spatial  location.  For  example,  data  layers  of 
vegetation,  rainfall,  soil  type,  evaporation,  and  surface-water  feaiures  can  be  used  to 
estimate  recharge.  These  spatially  referenced  data  layers  or  coverages  constitute  the 
input  for  distributed  hydrologic  models. 

Object  Oriented  Modeling  In  GIS 

Recent  advances  in  object-oriented  programming  have  begun  to  influence  the  field 
ofGISaswell.  The  recent  paradigm  isobjcct  orienlalion  in  GIS.  This  is  based  on  the 
same  promises  as  the  other  object  oriented  programming  languages.  Properiiessuch  as 
data  hiding,  encapsulation,  polymorphism,  and  Inheritance  ate  employed. 

In  object  oriented  programming,  data  and  procedures  are  treated  as  a single  entity 
referred  to  as  an '‘object."  Data  hiding  and  encapsulation  refer  to  the  ability  of  the 
process  to  hide  the  inlemal  structure  of  an  object  from  the  user.  With  encapsulation,  ills 
only  possible  to  manipulate  the  object's  data  using  a set  of  predefined  functions,  thus 
ensuring  data  independence.  The  internal  definition  of  the  data  structure  can  then 
change,  without  influencing  what  the  user  perceives. 


One  orihemoslimponinccoacepu  in  object-oriented  systems  Is  inheritance.  In 
object  oriented  programming,  inheritance  rerers  to  the  ability  to  create  a new  type  as  an 
extension  orthe  existing  type.  More  general  classes  Ihal  contain  the  generic  type  of 
object  am  brat  defined,  and  then  subclasses  are  created  as  specialiaalions  of  the  more 
generic  class.  These  subclasses  inherit  all  properties  ofthe  parent  class  and  add  some 
more  of  their  own.  For  instance,  when  emnting  an  object  class  to  represent  water 
network  nodes,  one  might  define  subclasses  to  be  a pump  and  reservoir,  which  would 
inherit  generic  water  node  characteristics  such  as  code  and  location,  along  with  their 
relationship  with  water  pipes,  and  aggregate  specific  characteristics  such  as  pump 
throughput  and  reservoir  capacity. 

Polymorphism  Is  the  other  impoitanl  feature  of  object  orientation.  This  refers  to  a 
program's  ability  to  use  objects  belonging  to  din'ereni  classes  In  a transparent  way, 
inlerpreiing  their  characteristics  on  the  run.  For  instance,  if  a program  needs  to  have 
access  to  an  attribute  called  code,  which  has  been  defined  for  several  diffetem  object 
classes  or  for  classes  with  common  ancestry,  it  is  not  necessary  to  know  the  object’s  class 
beforehand. 

The  advantages  of  this  recent  modeling  approach  are  numerous.  The  concept  of 
data  and  class  makes  it  possible  to  minimize  the  data  allocation  problem  by  categorizing 
data  into  classes.  The  property  of  inheritance  ensures  that  one  class  denotes  a number  of 
subclasses,  which  inherit  similar  properties  from  the  main  class.  Data  hiding  and 
encapsulation  make  it  possible  to  treat  dais  and  process  as  a unit,  thus  making  the  process 
mom  cflicient.  Polymorphism  makes  it  possible  to  create  different  data  features  that 
share  the  same  name,  thus  minimizing  the  memory  access  time  ofthe  process.  With  the 
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(ools  of  object  oriented  prognunming  and  advanced  geo*ptocessing  capabilities  of  CIS 
software,  the  integration  of  CIS  and  hydrologic  modelbig  enhance  the  modeling  pmcess 
significanUy. 

AppUcation  of  Optimization  Methods  in  Water  Resources 
The  management  of  water  resources  requires  an  integrated  approach  based  on  a 
careful  investigation  ofboth  the  hydrologic  factore  and  human  water  use.  Depending  on 
the  type  of  optimization  problem  to  be  addressed,  the  factors  to  be  considered  include  the 
type  of  engineering  adopted  in  the  design,  the  management  of  the  system,  economical 
and  legal  issues  concerning  water  rights  and  water  use,  ecological  factors  affected  by  the 
water  use  and  development,  and  social  factors  involving  the  decision  mokers  and 
stakeholders.  A combination  of  the  factors  considered  is  employed  to  formulate  a policy 
that  enables  managemenl  of  water  resources  on  a sustainable  basis. 

One  important  tool  used  in  formulation  of  decision  support  systems,  policy 
frameworks,  and  managemenl  plans  of  dtfferenl  scales  is  optimization.  There  are  a 
number  of  mathematical  tools  used  to  optimize  the  ulilization  of  water  resources  subjeci 
to  various  constraints.  The  constraints  in  most  cases  require  meeting  some  ecological 
needs,  quantitative  and  qualilative  hydrologic  factors,  and  a forecasted  demand. 

Optimization  is  a very  rigorous  process  that  requires  a very  careful  outlining  of  the 
following: 

• Dcfuiition  oflhe  problem. 

• Identification  of  paramelors. 

• Fomiulalion  of  the  objective  funclion. 

• Formulation  of  the  constraining  conditions. 

Definition  of  the  problem.  Problem  defimtion  requires  understanding  the  process 


that  needs  to  be  optimized.  At  this  stage,  a thoro 


; of  cause  and  efTeci  is 
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required  through  a detailed  analyaia  of  the  problem.  For  example,  a general  ^timiaation 
problem  in  water  resources  requires  idemification  of  the  water  sources,  investigation  of 
the  hydrological  and  consumptive  stresses,  and  various  options  that  could  bo  looked  into 
during  the  problem  formulation.  Detailed  historical  records  of  hydrologic  data  and 
consumption  records  need  to  be  analyzed,  and  forecasts  need  to  be  made.  Then  a 
decision  has  to  be  made  if  the  problem  is  a maximization  or  minimization  problem. 

Identincaiion  of  parameters.  The  parameters  are  the  basis  for  the  formulation  of 
the  objective  function.  The  patamelers  need  to  be  identilied  as  dependent  und 
independent  parameters. 

Formulation  of  the  objective  functlau.  Formulation  of  the  objective  function 
requires  writing  a series  of  linear  or  non-linear  mathematical  relationships  that  dellne  the 
pumping  requirements  subject  to  the  constraints  as  explained  earlier.  Choosing  the  most 
apprapriale  formulation  scheme  impacts  the  final  analysis  very  significantly.  Depending 
on  how  the  objective  function  and  the  constraining  functions  are  formutaied,  the 
optimization  problem  could  be  classified  as  integer  programming  (IP),  quadratic 
programming  (QP).  linear  programming  (LP).  non-linear  programming  (NLP),  mixed 
integer  programming  (MIP),  dynamic  programming  (DP),  and  so  on.  The  solution 
techniques  also  difTer  from  graphical  methods  used  to  solve  systems  of  linear  equations, 
non-linear  equations,  and  differential  equations.  Definition  of  the  solution  space  over 
which  the  solution  Is  feasible  is  important  in  some  cases.  This  region  is  generally 
bounded  by  the  constraining  equalionsoflhe  objective  function. 

FnrmulatioD  of  the  constrainlne  coodltloos.  The  eonatraining  conditions  dellne 
the  maximum  value  or  the  minimum  value  that  could  reasonably  be  achieved  without 
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jeopardizing  the  system.  For  example,  in  a water  supply  system  from  a river,  the 
maximum  amount  of  water  that  could  be  withdrawn  irom  the  river  could  be  constrained 
by  the  minimum  (low  requirement  that  needs  to  be  maintained  at  all  lime. 

DeCinltion  of  conrergeoce  criteria.  This  criterion  is  a mathematical  criterion  that 
defines  the  conditions  for  the  optimal  value  of  the  parameter  to  be  optimized  on  the  basis 
of  some  sets  oT mathematical  rules.  For  example,  ifhead  is  to  be  optimized,  then  the 
difference  between  the  head  value  at  one  lime  step  is  compared  with  the  head  at  the  next 
lime  step  and  the  difference  compared  with  the  convergence  criterion.  This  is  applicable 
in  a more  complicated  optimization  scenario. 

Optimization  modeling  is  a general  class  of  problems,  which  require  minimization 
or  maximization  subject  to  a set  of  coitstroints.  A number  of  optimization  models  have 
been  developed  to  address  the  issue  of  groundwater  management.  Most  models  are 
developed  to  formulate  optimal  pumping  strategies  by  employing  cost  minimization  or 
pumping  maximization  subject  to  water  quality  and  quantity  constraints 

Optimization  algorithms  developed  for  management  of  groundwater  systems  have 
to  be  coupled  with  a simulation  model.  The  coupling  of  single-objective,  linear 
groundwater  optimization  models  is  broadly  classified  on  the  basis  of  the  approaches 
used.  The  two  approaches  are  embedded  and  response  matrix  approaches  (Gorelick. 
1983). 

Embedding  methods  involve  development  of  a linear  set  of  equations  known  os 
response  equations  from  the  partial  diRercnIial  equations  of  groundwater  flow  and 
Irunsport.  The  solution  of  the  resulting  equations  is  obtained  by  using  numerical  methods 
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such  as  llnile  difTercncc  or  finite  elements.  Applications  of  this  approach  have  been 
successfully  tested  by  Aguodo  and  Remson  (1974)  and  Willis  and  Liu(l9S4). 

The  response  matrix  approach  is  developed  fiom  hydraulic  or  watcrK|uality  models 
ofthe  aquifer  system.  The  simulation  model  Is  used  to  generate  the  response  of  the 
aquifer  systemic  pumping  or  recharge.  The  resulting  matrix  of  coefficients  is  used  to 
convert  pumping  stresses  into  linearized  changes  in  drawdovms  and/or  changes  in 
concentrotion  ofthe  ambient  water.  The  set  of  unit  responses  is  formulated  for  individuaJ 
wells  and  superimposed  to  model  multiple  wells,  which  foim  a set  of  linear  equations  that 
is  solved  systematically.  Hovrever,  if  the  problem  is  formulated  as  a non'linear  problem, 
then  the  solution  is  obtained  iteratively.  The  response  matrix  approach  has  been  used  to 
solve  multi-objective  optimization  problems,  (e.g.,  Louie  el  al..  1 984  and  Yeh  ei  al.. 
1992). 

Research  Objectives  and  Methodology 

Them  is  a considerable  amount  of  knowledge  that  describes  the  process  of  saltwater 
intrusion.  In  the  past  couple  of  decades,  different  mathematical  approaches  hsve  been 
developed  to  model  the  pioccss  of  seawaler  intrusion.  There  are  also  well-documented 
research  results  and  a number  ofGCMs  that  make  climatic  variability  forecasts,  usually 
on  the  time  scale  of  a century. 

The  majority  of  research  on  saltwater  intrusion  is  primarily  focused  on  assessing 
the  effect  of  freshwater  withdrawals  from  coastal  aquifers.  Although  it  has  been 
established  that  climate  change  impacts  influence  Ihe  process  of  seawater  intrusion,  there 
is  0 gap  in  the  dynamic  Unking  of  climate  change  impacts  to  seawater  mtnision  and 
quantification  of  the  flow  and  transport  processes.  There  Is  also  a gap  in  the  development 
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of  on  integrated  approach  to  manage  coastal  water  resources  by  incorporating  the  impoels 
of  climate  change  and  fipshwater  withdrawal. 

The  objectives  of  the  research  described  in  this  dissertation  were: 

• To  obtain  a better  understanding  of  the  impacts  that  climate  change  and  sea  level 
rise  will  have  on  seawater  Imnision  In  coastal  aquifers. 

• To  bridge  the  gap  ofknowledgc  in  the  integrated  assessment  of  coastal  aquifers, 
particularly  management  on  a lime  scale  of  a century. 

« To  investigate  the  sensitivity  of  coastal  aquifer  parameters  and  hydrologic  stresses 
to  Che  process  of  saltwater  intrusion  coupled  with  the  impacts  of  sea  level  rise. 

The  research  objectives  described  above  address  a wide  array  of  subtopics.  As  a 
result,  each  objective  required  Ihe  application  of  a combination  of  hydrologic, 
hydrogcologic,  and  spatially  distributed  numerical  modeling  methods. 

Investigation  of  Ihe  impacts  of  sea  level  rise  involved  transient  modeling  overs 
given  area,  and  as  a result,  a CIS  approach  was  chosen  to  perform  the  task.  Investigation 
of  the  impacts  of  sea  level  rise  and  climate  change  on  Ihe  groundwater  fknv  regime  io 
coastal  aquifers  required  development  of  a groundwater  flow  model  that  incorporates  the 
transient  nature  of  sea  level  rise  as  a boundary  condition.  Similarly,  investigation  of  the 
management  of  freshwater  withdrawals  required  application  of  optimization  methods  to 
develop  optimal  piunpmg  strategies.  Also,  bridging  the  gap  of  knowledge  In  Ihe 
Imegraicd  management  of  coastal  aquifers  required  development  of  a procedure  that 
combines  sea  level  rise  modeling,  groundwater  flow  slmulabon.  and  optimizalion. 

In  order  to  proceed  with  the  research,  three  major  methodologies  were  identified, 
viz.,  CIS  for  sea  level  rise  modeling,  numericsl  modeling  for  groundwater  How 
simulation,  and  optimization  for  management  of  freshwater  withdrawals.  Afier 


idenliryinglhe  meihodologies.  detailed  procedures  tor  each  methodology  were  outlined. 
The  general  classes  of  methodologies  used  in  this  research  are: 

■ Development  ota  CIS  based  sea  level  rise  modeling  method. 

• Development  of  a variable  density  groundwater  flow  modeling  computer  code  that 
incorporates  the  effects  of  sea  level  rise  by  modifying  an  esisting  three- 
dimensional,  variable  density,  groundwater  flow  simulation  code. 

• Development  of  an  optimization  model  based  on  evolutionary  algorithms. 

• Development  of  an  approach  for  integrated  coastal  water  resources  management  by 
combining  sea  level  rise  modeling,  groundwater  now  simulation,  and  optjmizaiion 
in  one  shell  program. 

As  pariofthis  research,  a number  of  computer  codes  were  developed  for  each 
methodology.  A shell  program  that  incorporates  all  the  three  components  also  was 
developed.  The  programming  languages  used  in  the  development  of  the  computer  codes 
include  Avenue.  FORTRAN,  and  Visual  Basic. 

Sea  level  rise  modeling  was  performed  for  coastal  areas  of  south  Florida. 
Groundwater  flow  simulation  was  peTformed  for  a coastal  aquifer  in  south  Florida  and  a 
coastal  aquifer  in  Turkey.  Optimization  modeling  was  performed  fora  test  case  ota 
coastal  aquifer  under  salinization. 

Dissertation  Organization 

The  first  chapter  of  this  dissertation  describes  the  general  background  required  for 
the  simulation  and  optimization  of  water  resources  management  problems  and  the 
application  of  GIS  in  water  resources.  The  attempt  is  lo  give  the  reader  a basic 
undersianding  of  the  topics  rather  than  a detailed  description  of  methods  and  concepts. 

The  process  of  seawater  imrusion  in  coastal  aquifers  is  described  in  the  second 
chapter  with  special  emphasis  on  numerical  modeling.  The  general  mathematical 
formulation  ofa  voriablo  denaity  groundwater  flow  problem  and  the  numerical  solution 
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ofihe  problem  are  discussed.  The  cnecis  lhal  climste  change  and  sea  level  rise  have  on 
ihc  groundwater  flow  rcgiine  are  also  brieHy  discussed.  Also,  ihe  numerical  three- 
dimensional  groundwater  flow  model  SEA  WAT  is  described. 

Chapter  three  describes  optimiaalion  methods  based  on  evolutionary  compulation. 
Specifically,  Genetic  Algorithms  (CAs)  ore  discussed,  and  their  application  and 
suitabiliQ^  in  search  and  optimization  ace  investigated. 

In  chapter  four,  a GfS  model  is  developed  to  investigate  the  impacts  that  sea  level 
rise  has  on  coastal  areas.  The  model  framework  is  described  in  detoil,  and  the  computer 
codcdevelopmentprocess  is  explained.  In  addition,  the  IPCC  scenarios  are  analyzed 
with  special  emphasis  on  global  ses  level  rise  and  on  methods  to  estimate  relative  sea 
level  rise  ofan  area  from  global  sea  level  rise  estimates. 

Chapler  live  describes  Ihe  method  adopted  to  incoeporafe  the  impacts  of  sea  level 
rise  In  the  modeling  of  groundwater  Ilow  in  coastal  aquifers.  The  boundary  conditions 
commonly  used  in  groundwater  flow  modeling  are  discussed,  and  modifications  to  a 
three-dimensional  variable  density  groundwater  flow  model  arc  presented. 

An  optimizalion  model  based  on  genetic  algorithms  is  developed  and  described  in 
chapter  sis.  Model  components,  solution  techniques,  and  state  and  decision  variables  are 
discussed  in  detail.  The  coding  aspects  of  the  new  groundwater  flow  optimizalion  model 
QENOPT,  which  Is  based  on  genetic  algorithms,  also  are  explained. 

Results  ofihc  simulation  and  optimization  modeling  are  presented  In  chapler  seven. 
The  shift  of  Ihc  saltwater-freshwater  interface,  simulated  heads  and  concentration  values, 
and  optimized  model  parameters  are  used  to  explain  the  results. 


Chapler  eight  is  the  final  chapter,  where  the  conclusions  drawn  Erom  the  research 
are  summarized  and  presented.  Suggestions  and  recommendations  for  future  research 
opportunities  also  are  presented. 

Appendix  A presents  the  concept  ofintegraied  modeling  as  applied  to  the 
management  ofcoasial  aquifers  subject  to  the  impacts  of  climate  change.  The  shell 
program  that  incorporates  the  three  modules,  viz.,  GIS.  simulation,  and  optimization,  is 
presented.  A discussion  ofsofiware  anddata  requirements  and  a series  of  transient  plots 


also  are  included. 


CHAPTER  2 

SEAWATER  INTRUSION  IN  COASTAL  AQUIFERS 
Concepts  and  Methods 

ScawQler  intrusion,  also  known  as  saliwoier  intrusion,  is  described  as  the  mass 
transport  of  saline  waters  into  zones  of  an  aquifer  that  were  previously  occupied  by 
freshwater  (Stewart,  1 999).  Degradation  of  groundwater  quality  due  to  seawater 
intrusion  is  characterized  by  increased  levels  of  salinity  that  exceed  acceptable  values  for 
potable  and/or  agricultural  purposes.  Seawater  intrusion  primarily  aifects  coastal 
aquifers,  which  serve  as  major  sources  of  freshwater  supply  in  many  parts  of  the  world. 
Estimates  show  that  coastal  plains  are  inhabited  by  about  7(«4  oflhe  world’s  population. 
The  steadily  increasing  demand  for  potable  water  and  the  threat  posed  by  saltwater 
intrusion  have  become  serious  issues  of  concern.  As  a result,  management  of  freshwater 
resources  in  coastal  aquifers  has  emerged  as  a challenging  problem  to  planners,  policy 


The  rate  of  mass  transport  of  seawater  and  its  pothway  in  the  subsurface  is 
determined  by  the  hydraulic  gradient  and  distribution  of  hydraulic  conductivity  (Stewart, 
1999).  The  flow  of  saline  water  into  the  frvshwater  aquifer  and  the  discharge  of 
freshwater  into  the  coastal  boundary  is  a long-term  process  that  takas  a considerably  long 
lime  to  reach  a new  state  of  equilibrium  (Van  Dam,  1 999).  The  dynamics  of  flow  in  the 
aquifer  are  governed  by  the  existence  of  pressure  gradients,  which  in  turn  are  determined 
by  aquifer  boundary  conditions  and  soumes  and  sinks.  Under  natural  ciicumslances.  the 
boundary  conditions  arc  determined  by  components  of  the  hydrologic  cycle  and  relative 
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sea  level  rise.  However,  anlhropogenic  faciors  such  as  increased  pumping  allcr  Uie 
process,  and  lead  to  seawoier  inlnisiort 

Salinity  of  walcr  Is  usually  measured  in  terms  of  either  the  Touil  Dissolved  Solids 
(TDS)  or  the  concemration  of  chloride  ions  (Cf).  Total  dissolved  soiids  (TDS)  describes 


principal  constituents  are  usually  calcium,  magnesium,  sodium,  and  potossiuni  cations 
and  carbonate,  hydrogencarbonalc  (HCO'sX  chloride,  sulfate,  and  nitrate  anions  (WHO, 

1 9961.  Chlorides  are  distributed  widely  in  nature  os  .salts  of  sodium  (/VoCT).  potassium 
(AlC/),  and  calcium  (CoC/;).  Tables  2-1  and  2-2  present  a classification  based  on  both 
methods.  The  data  presented  in  Table  2-1  was  obmined  from  Chemical  Rubber  Company 
(CRC,  1982) 


Table  2-1.  Properties  of  various  concentrai  onsof  seawater. 


Salinity 
(parts  per 
thousand) 


Chlorinity 
(pans  per 
thousand) 


Absolute 
viscosity 
(centipoises,  cP 
at  20  Tl 


Kinematic 
viscosity 
(centistokes,  cS 
ai20°g 


Tbe  Problem  of  Seawater  Intrusion 


The  problem  of  seawater  intrusion  is  characterized  by  the  presence  of  two  fluids  of 
di/ferent  solute  concentration,  freshwater  and  saltwater.  There  are  two  difierant 
approaches  to  the  solution  of  the  probiein  of  seawater  intrusion.  The  first  approach 
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kncpwn  as  chc  sharp  interface  model  involves  assmning  the  two  fluids  are  immiscible  and 
that  there  exists  a sharp  interface  that  separates  them.  The  second  approach,  known  as 
the  variable  density  model,  considers  the  two  fluids  as  miscible  and  incorporates  a 
transition  loae  that  separates  the  freshwater  zone  with  the  saltwater  zone.  This  approach 
calls  for  a method  that  accounts  for  variability  of  density  in  the  flow  zone. 

In  a situation  where  the  interface  is  approximated  to  be  sharp,  and  with  a 
simplifying  assumption  of  a predominantly  horizontai  flow,  the  problem  of  finding  the 
interface  may  sometimes  be  solved  Inclosed  foim'.  However,  neglecting  the  existence 
of  hydrodynamic  dispersion  and  a mixing  zone  is  generally  inadequale  in  terms  of 
accurately  representing  the  flow  regime.  As  a result,  most  of  the  recent  two  and  three- 
dimensional  models  of  seawater  intrusion  account  for  the  transilion  zone  and  vertical 


The  mathematical  formulation  of  the  coupled  transport  and  flow  terns  that 
represent  the  governing  partial  difierenlial  equations  in  three  dimensions  is  generally 
solved  using  numerical  methods.  The  major  diffieuity  involved  in  finding  analytical 
solutions  is  the  non-linearity  of  the  parameters.  There  are  some  quasi-analylical  and 
analytical  solutions  for  one  and  two-dimensional  problems. 

The  numerical  melhods  used  in  solving  seawater  intrusion  problems  gcneraJly  use 
either  finite -difference  or  finite-element  methods  and  solve  for  hydraulic  head  and  solute 
concemraUon.  The  hydraulic  head  is  customarily  known  as  environmental  head  as  it  can 

' An^iulion  is  said  la  be  s close^fonn  solution  if  it  sotvesagiven  pmblefn  in  terms  onuncttoiis  and 

not  be  consideted  closed-form.  However,  the  choice  of  what  to  call  closed- form  and  wtot  not  Is  lalta 
aibilrary  since  a new  "closed-fbmi''  fiincllon  enuld  simply  he  defined  in  terms  of  the  infinite  sum. 


be  head  either  in  the  Iresh  or  saliwalef  zoneofihc  aquifer  and  can  readily  be  measured  by 
using  a piezometer.  Many  solution  lechniques  solve  for  the  hydraulic  head  in  terms  of 
the  equivalent  freshwater  head,  which  can  then  be  converted  to  environmental  head  by 
using  the  concentration  and  density  terms  as  required.  The  equivalent  freshwater  head  is 
a measure  of  the  environmental  head  in  the  sahwater  zone  representing  the  level  to  which 
the  water  level  w ould  rise  if  the  piezometer  were  filled  with  frcdiwaler. 

The  primary  concern  to  studying  the  problem  of  seawater  intrusion  is  the 
permanent  loss  of  a freshwater  aquifer  due  to  excessive  levels  of  salt  in  the  aquifer.  Once 
an  aquifer  is  subject  to  intrusion,  reversing  the  process  and  reclaiming  the  aquifer  is  a 
very  difricult  and  expensive  undertaking.  There  have  been  some  attempts  to  mitigate  the 
impacts  of  aquifer  loss  to  seawater  intrusion  in  the  western  coast  of  the  United  Slates. 
Two  such  examples  arc  Ihe  projects  in  Orange  and  Los  Angelos  counties  In  California, 
where  there  are  efforts  underway  to  contain  the  advance  of  seawater  intrusion  by  creating 
subsurface  hydraulic  barriers. 

Malbematical  Modeling  of  Seawater  Intrusioo 

The  pioneering  work  in  the  investigolion  of  saltwater  intrusioo  was  done  by  Badon 
Ghyben  (1888)  and  Hertzberg  (1901).  The  Ghyben-Hcrtzbcrg  formula,  as  it  has  become 
to  be  known,  established  a relationship  between  elevation  of  water  table  in  an  unconfined 
aquiferlo  the  elevation  of  the  boundary  known  as  soltwater-freshwalcr  interface. 

Looking  at  Figure  2-1,  Ihe  relationship  between  positron  of  the  interface  and  the  wafer 
table  head  is  written  os: 


(2-1) 


where  p,and  ^.aie  densities  of  fresh  and  salt  water  respectively.?  is  the  depth  of  the 


interface  below  Ihe  measuring  point,  normally  mean  sea  level,  and  h is  the  elevolion  of 
the  water  table.  Assuming  a density  of  1.025  kg/m’  for  the  saltwater  and  1,000  kg/m’  for 
the  freshwater,  the  depth  of  the  interface  is  approsimsiely  40  times  the  value  of  h.  This 
serves  as  a preliminary  estimate  for  the  position  of  the  interface  in  many  cases. 


Sharp  Interface  Approach 

The  simplified  assumption  of  a sharp  inlerface  is  by  no  means  the  most  accurate 
lepresemation.  When  Ihe  thickness  of  the  freshwater-saltwater  transition  zone  is  small 
compared  to  the  thickness  of  the  aguifer.  then  it  Is  possible  to  assume  that  the  treshwaler 
and  saltwater  domains  are  separated  by  a sharp  interface  (Es.said.  1009).  Thecwo 
domains  of  saltwater  and  freshwater  are  coupled  by  Ihe  inierlacial  boundary  condition  of 
continuity  of  flux  and  pressure  (Essaid.  1999).  In  situations  where  the  thickness  of  the 
Inlerface  is  significant,  the  shaip  interface  ceases  to  exist,  leading  to  a mixing  zone  where 
hydrodynamic  dispersion  lakes  place.  Although  a sharp  interne  approach  reproduces 
Ihe  legionaJ  flow  dynamics  end  the  response  of  Ihe  sharp  intcriuce  to  stresses,  it  does  not 
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provide  informalion  on  thenaiuieofilie  iransilion  zone.  In  three-dimensions,  the 
inlerfaciol  bountUty  condition  is  hii-hly  non-linear  (Bear,  1979),  as  a resnlt.  both  the  How 
and  the  aquifer  have  to  be  aasumed  to  be  horizontal  in  order  to  find  a solution  (Essaid, 

1909).  In  orderlo  obtain  the  ihree-ditnensionaJ  solution,  the  horizontal  flow  is  integrated 

vertically, 

Huhben  (1940)  developed  equations  that  related  densities  of  freshwater  and 
saltwater  with  the  elevation  of  a sharp  interface  in  tenns  freshwater  heads  using  the 
concept  of  a potential  function.  The  freshwater-saltwater  relationship  was  developed  by 
equaling  pressure  at  a point  on  the  inlerface.  Equaling  the  pressure  at  the  interface  for 
both  the  freshwater  and  saltwater  eovironments,  the  freshwater  head  is  computed  as: 


and  the  saltwater  head  is  computed  as: 


where  h/ond  A,  are  freshwater  and  saltwater  heads,  p^and  p,are  densities  offresh  and 
sail  water  re^eclively.  and  : is  the  elevation  of  the  inlerface  above  a datum.  The 
schematic  description  of  the  approach  formulated  by  Hubbcrt(l940)and  modified  by 
Reilly  and  Goodman  (1985)  is  presented  in  Figure  2-2. 

Two  eases  are  considered  in  Figure  2-2  to  represent  the  motion  of  freshwater  and 
saltwater.  Both  the  freshwater  and  saltwater  are  assumed  to  be  two  Immiscible  fluids  of 
different  densities.  Case  I represents  the  situation  when  one  fluid  flows  while  the  other  is 
static,  while  Case  2 depicts  the  condition  when  both  fluids  arc  in  motion.  Notation  of  the 
subscripts  adopted  in  the  sketch  is  such  that  for  any  parameter  trj,  the ; represents  the 


fluid  cype  while;  represenls  the  parameter  number.  Thus,  ihij  is  the  potential  of 
freshwater  while  <I>2,  is  saltwater  potential.  Similarly,  qi  Is  freshwater  flow  vector,  while 
q2  is  saitwolcr  flow  vector.  The  angle  of  inclination  of  the  interface  is  represented  by  o 
and  : represenls  the  elevation.  For  the  system  to  be  in  equilibrium,  Ihe  saltwater  has  to  be 
stationary  {Case  I).  and  under  equilibrium  conditions,  the  elevation  of  a point  on  the 
interface  above  a datum  could  be  determined  by  Ihe  formula; 


Figure  2-2.  Description  of  the  potential  surface  approach. 

The  sharp  interface  approach  is  hirther  clas.siricd  into  two  groups  as  illusiraled  in 
Figure  2-2.  The  first  case  assumes  the  saltwater  to  be  static  while  the  freshwater  is  In 
motion.  This  essentially  reduces  the  problem  to  a one  fluid  How  with  a free  surface 
boundary'  condition.  The  saltwater  is  assumed  to  move  in  respoitse  to  the  flow  of 
freshwater.  Glover  (1 959)  developed  a sharp  interface  model  to  study  the  penem  of 
freshwater  flow  in  a coastal  aquifer.  Bennett  and  Oiusli  [1971)  investigated  coastal 
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BTOundwatcr  flow  in  Ponce,  Puerto  Rico,  uiing  the  saine  approach.  The  saltwater  beneath 
an  oceanic  island  was  studied  by  Fetter  (1972)  assuming  sharp  interface.  Similai  studies 
include  Voss  (1984a).  Guswa  and  Leblanc  (1985),  and  Reilly  et  ai.  (1987b). 

The  second  group  of  sharp  interface  problems  assumes  that  both  freshwater  and 
saltwater  ate  in  motion  as  described  in  Case  2 of  Figure  2-1  The  two  fluid  approach 
requires  the  flow  system  to  be  ueated  as  a two  separale  systems  with  a common  free- 
surface  boundary  (Reilly,  1993).  The  solution  ofthese  types  of  problems  requires 
solution  of  two  dilTercm  partial  differential  equations  that  govern  the  flow.  The  position 
ofihc  interface  is  determined  by  equation  2-4  (Hubbeit,  1940).  The  governing  sets  of 
mass  balance  equations  are  written  in  lernis  of  the  freshwater  (equation  2-5)  and  saltwater 
(equation  2-6)  are: 

Sk, 

= (2-5) 

„ Bh 

= (2-6) 

’-(ll*  (I),*  (I), 

directions,  his  head  (L).  q(LT')  is  specific  discharge.  S is  specific  storage  (L'').  g is  a 
sini/source  term  (T').  The  specific  discharge  for  both  freshwater  and  saltwater  is 
expressed  in  terms  ofDarey's  equation  of  How  for  constant  density  fluids  as: 

g,=-KVh,  (2-7) 


q,=-KVh, 


(2-8) 
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The  second  case  of  the  shaip  interface  approach  was  used  to  model  saltwater 

iiurusioninthcsouthempanofLongIsland,  New  Yoric(Pindcr  and  Page,  1977).  A 

finiie-difrerenee  solution  of  the  sharp  interface  with  two  flow  fields  was  also  developed 
10  model  area!  flow  of  freshwater  and  saltwater  separaled  by  an  inlerfacc  (Mercer. 
Larson,  and  Faiisl,  1980).  Similariy,  a sharp  interface  quasi-thiee-dimensional  numerical 
model  SHARP  was  developed  and  solved  using  Hnile.dif1crence  simulation  (Essaid, 

1 990a  and  1 990b). 

V ariable  Density  Approach 

The  variable  density  conceptualization  of  the  saltwaier-lreshwater  problem  is  based 
on  the  mixing  of  the  two  fluids,  thus  influencing  the  density  and  viscosi^  of  the  fluid. 
Due  to  the  variability  in  density  of  the  fluid  as  a result  of  diHerent  levels  of  concentration 
indifferent  zones,  the  flow  is  essentially  a variable  densi^  flow  problem.  Solution  of  a 
variable  density  flow  problem  normally  requires  solving  the  non-linear  mass  balance 
equations  for  fluid  and  salt.  The  governing  equations  are  solved  for  fluid  pressure  or 
hydraulic  head  and  concenuntion.  In  order  to  find  the  solution,  the  fluid  is  assumed  to  be 
incompressible.  Isothermal  conditions  also  are  assumed  for  the  majority  of  flow 
problems  in  varying  hydrogeologic  settings. 

The  variable  density  system  conceptualization  used  in  ihe  development  of  the 
numerical  model  SUTRA  (Voss,  1984b)  is  presented  to  describe  a variable  density  flow 
formulation.  The  primary  step  in  the  formulation  of  a variable  density  model  is 

conservation  of  mass.  The  fluid  mas.s  balance  is  expressed  as  the  sura  of  pure  water  mass 

and  solid  malrix  mass  as: 


^^  = -V(«pv)+S 


(2-9) 
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where  II  is  porosily  (dimensionless),  p is  fluid  densily  (ML‘‘),  0 is  fluid  mass  soumc 
(ML  ’r'),  V is  average  fluid  velocity  (LT')obtained  by  dividing  the  speciflu:  dischnige? 
by  the  porosity  n.  and  l is  lime  (T).  The  average  fluid  velocity  depends  on  fluid  pressure 
and  density.  Using  Darcy's  law  and  variable  densily.  it  is  expressed  as: 

► = (2-10) 
where  k is  the  aquifer  permeability  (L^,  p is  the  fluid  viscosity  (ML‘'T'),  and  gisthc 


g = -lg[7(e/evo«bnJ  (2-H) 

In  order  to  make  use  of  the  fluid  potential,  the  hydraulic  conductivity  AT  is  defined 
for  a constant  density  fluid  as: 


(2-12) 


Substituting  the  constant  density  Darcy  equations  2-6  and  2-7  into  equation  2-8  yields 

^ (VA'-Arg)J+0  (2-13) 

The  left  hand  aide  of  the  mass  balance  equation  can  also  be  expressed  in  terms  of 
pressure  and  concentration.  But  first,  the  storativiiy  of  a saturated  aquifer  is  related  to 
pressure  and  defined  as: 


BP 


(2-14) 


where  the  term  = (I -n)or+n^  is  the  specific  pressure  slorativity  (LT^M  ’),  o istl 
aquifer  matrix  compressibility  (LT^W'),  and  gisthc  fluid  compressibility  (LT^M"'). 


y~y  a a 

.equalion2-ll  yields: 


(2.17) 


(2-20) 

(2-21) 


D„  = D„  =!(</, -rf,X-.v.)  (2-22) 


■(LT’); 


>V  = magnitude  of  vertical  componenl  of  v(LT'); 
di  = the  longiludinal  dispersion  coefncicnt  (L'T');  and 
dr  =■  the  transverse  dispersion  cocfncienl  (LV). 


(2-23) 


(2-24) 


oi  = ihc  longiludinai  dispersivlty  of  the  solid  mairix(L):  &nd 
ffr*  die  transverse  dispersivity  of  the  solid  malrix  (L). 

A three-dimensional  version  of  the  variable  density  groundwater  flow  model 
SUTRA  was  released  in  2003  (Voss  and  Provost,  2003).  The  new  code  still  maimains 
the  basic  features  of  the  two-dimensional  model  (Voss,  1984)  and  simulates  fluid 
movement  and  the  transport  of  either  energy  or  dissolved  substances  in  a subsurface 
environment.  The  new  code  employs  a two-  or  three-dimensional  finiie-eiement  and 
finite-difference  method  to  approximate  the  governing  equations  that  describe  the  two 
imerdcpendempracesses  that  are  simulated.  The  simulated  processes  am  fluid  density- 
dependent  saturated  or  unsaturated  groundwater  flow  and  solute  transport  or  transport  of 
thermal  energy  in  groundwater  and  solid  matrix  oflhe  aquifer. 

Huyakom  ei  al.,  (1987)  developed  a three-dimensional  ^niie  element  model  to 
simulalc  saltwaler  inlrusion.  The  governing  equations  of  the  model  consisted  of  a 
variable  den^ty  fluid  flow  equation  and  an  advcctive-dispcrsive  transport  equation.  The 
frcshwaier  hydrruilic  head  and  concentration  are  treated  as  dependent  variables.  The 
three-dimensional  miscible  flow  model  is  written  cut: 


where  = the  hydraulic  conductivity  tensor  (LT');  A = reference  hydraulic  head 
referred  to  as  the  freshwater  head  (L);  x,  (i=l,2,3)  = Canesian  coordinates  (L);  7 
density  coupling  coefficient;  C = solute  concentration  [flS  C S l)(ML  ');  e,  =/" 


component  of  the  graviialional  1 


l.ej-0).  ;Ss  specified  storage 
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[L  );  T = lirne  (7);^  = porosily;  q = voluinetric  flownuc  of  the  sources  per  unit  volume 
ofthe  porous  medium  p and  density  of  the  mixed  fluid  (freshwater  end 

seawater)  and  reference  (freshwater) density  respectively  (MV‘).  Itinust  be  noted  that 
the  solute  concentration  is  zero  at  the  minimum  (reference)  density. 

The  reference  bead  and  dertsity^coupling  coefficient  in  Equation  2*1  are  defined 


(2-26) 


n = ~-  (2-27) 

where  p = fluid  pressure  « gravitational  acceleration  (Wir');  2 * elevadon 

above  datum  iL);  e = density  diflercncc  ratio  defined  as: 


and  C,  = solute  concentration  [ML  ’)  corresponding  to  the  maximum  density  p,  (ML-’). 
The  actual  hydraulic  conductivity  in  (1)  which  depends  on  fluid  density  and  viscosity, 
is  generally  a function  of  solute  concentration.  It  is  defined  as: 

(2-2fl) 

where  i,' intrinsic  pcimonbilily  tensor  (l');  /r » dynamic  viscosity  of  the  fluid 
(S'  = S'j)  (Afi‘'r');and  p,  * viscosi^  of  the  freshwater  (Mt'r'). 

ITte  density  of  the  fluid  in  the  mixed  region  is  treated  as  n linear  fiinaion  of 


concentration  and  can  be  expressed  os: 


(2-30) 


The  odvKiive-dispersive  equalion  used  lo  describe  Ihe  sail  ironsfcn  is  expressed  in 
the  following  form: 

o->i> 


rvhere  Z)^=^D^.wilh  Dt  being  Ihe  dispersion  tensor  (2,^r')defined  by  (Bear.  1979),  V, 
" Darcy  velocity  vector  {IT' ). 

The  Darcy  velocity  vector  is  expressed  as: 


-u/Ce,  I 


C-32) 


where  Kl  - the  hydraulic  conductivity  at  ibe  reference  (freewater)  condition  {LT'). 

Henry  (1964)  developed  a semi-analytreal  model  that  dedoesthe  location  and 
shape  of  the  saltwaler-freshwaier  interface  under  the  condition  of  a constant  seaward  flux 
of  freshwater  toward  an  oceanic  boundary.  This  semi-analytical  model,  in  which 
dispersion  is  assumed  to  be  constant,  still  scrve.s  asone  ofthe  standard  benchmark 
problems  for  testing  newly  developed  numerical  models. 

There  has  been  a significant  improvement  in  the  modeling  of  variable  density 
groundwater  flow.  This  development  was  greatly  enhanced  by  the  development  of 
numerical  models  and  Ihe  advent  of  computers.  The  numerical  models  currently  in  use 
include  the  U5.  Geological  Survey  two-dimensional  finite  element  model  SUTRA 
(VOSS,  1984b).  the  three-dimensional  finite  difference  model  HST3D  (Kipp,  1997), 
MOCDENSE  (Sanford  and  Konikow.  1985),  and  SEA  WAT  (Guo  and  Langevin,  2002). 
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In  letms  of  Ihe  representation  of  the  saltwalcc-freshwaier  imccfacc,  the  generaliyed 

(Bear.  1 999).  These  are: 

I . Three-dimensional  sharp  interface  models: 

SEAWAT:  Variable  Density  Groundwater  Flow  Model 

density,  finite  difference  model  called  SEAWAT  (Guo  and  Langevin,  2002).  SEAWAT 

MODFLOW  (McDonald  and  Haibaugh.  1988)  and  MT3DMS  (Zheng  and  Wang.  1998). 
Similar  to  MT3DMS,  SEAWAThas  Ihe  capability  to  use  a variety  of  solution  techniques, 
and  it  maintains  the  modular  structure  of  MODFLOW  and  MT3DMS.  As  a result.  Input 
files  developed  for  MODFLOW  and  MT3DMS  can  be  used  directly  as  input  files  for 
SEAWAT  with  only  lilUe  modification.  A version  of  SEAWAT  based  on  MODFLOW 
2000  was  released  very  recently  (Langevin.  Shoemaker  and  Quo.  2003).  SEAWAT-2000 
was  designed  by  combining  a modified  version  of  MODFLOW-2000and  MT3DMS  into 
a single  computer  program.  The  concept  of  process  used  in  MODFLOW-2000  was  also 
used  b Ihe  development  of  SEAWAT-2000.  A process  is  defined  as  pan  of  Ihe  code  that 
solves  a fundamental  equation  by  a specified  numerical  method.  SEAWAT-2000 
contains  all  of  Ihe  processes  distributed  with  MODFLOW-2000  and  also  includes  the 
variable-density  flow  process  and  the  blegraled  MT3DMS  transport  process. 

The  development  of  SEAWAT  required  modifying  MODFLOW  so  that  it  is  based 
on  Ihe  principles  of  mass  conservation.  Instead  of  Ihe  original  formulation  of  fluid 


I of  using  environmental  heads.  SEAWAT  i 
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equivalent  fteshwaier  heads  (Lusczynski,  1961).  In  order  to  understand  the  concept  of 

equivalent  freshwater  head,  consider  two  piezometers  as  shown  in  Figure  2-3.  one  filled 

with  freshwater  and  the  other  filled  with  saltwater.  The  height  of  the  water  level  in 

piezometer  I above  the  datum  is  where  p,  is  the  density  of  Ireshwatet.  The 
P,t 

freshwater  bead  at  the  measuring  point  (bottom  of  the  piezometers)  is  equivalent  to  the 
elevation  of  the  water  level  in  piezometer  1 . which  is  expressed  using  the  energy  equation 
as  described  in  Equation  2-2.  Similarly  the  head  in  the  saltwater  piezometer  is  expmssed 
as  Equation  2-3.  Where  Ay  is  equivalent  freshwater  head  (L),  h is  head  (L),  fls  pressure 
(ML‘'T’),  p,  is  freshwater  density  (ML’’).  is  saltwater  density  (ML'“).  g is 
acceleration  due  to  gravity  (LT^),  and  Z is  elevation  above  datura  (L). 

The  basic  assumptions  involved  in  the  development  of  SEA  WAT  are  based  on 
Darcy’s  law  for  flow  through  porous  media.  Pick's  law  of  diffusion  for  the  dispersive 
Transport,  the  prevalence  of  isothermal  conditions,  and  applicabilirv  of  the  specific 
storage  coneept  to  confined  aquifers.  The  porous  medium  is  assumed  to  be  fully  saturated 
with  water,  and  a single,  fully  miscible  liquid  phase  of  very  small  compressibility  is  also 


assumed  to  fill  the  porous  medium. 


Piezcrrvter  1 - Freewater 


Figure  2-3.  Schemalic  repmenuiion  of  the  concept  of  equivalent  freshwater  head 
The  process  of  diffuaion  enables  particles  to  be  transferred  Irotn  zones  of  hl^ 
concentration  to  zones  of  low  concentration.  Pick's  law  of  diffusion  states  that  the  mass 
flux  of  particles  is  proportionate  to  the  concentration  gradient; 

«=-d,7C  (2-33) 

where  ^ is  the  mass  flux  of  particles  and  tfo  Is  the  molecular  diffusion  coefficient  in 
water.  The  classical  dispersion  theory  has  been  developed  by  Tnylor  (1933).  Josselln  de 
Jong  (1 958).  Saflinan(1959,l960),Scheiidcgger  (1960).  Bear  and  Bachmat  (1967).  and 
Fried  and  Combamous  (1971).  The  mathematical  representation  of  the  dispersive 
transport  process  Is  analogous  to  that  of  Ficlc’s  law  and  accounts  for  the  mixing  Involved. 
With  the  dispersion  coefficient  of  the  medium  defined  as  D.  the  mass  flux  is  expressed 


<t)  = -£l7C 


(2-34) 


The  head  values  calcuiaicd  by  SEA  WAT  are  equivalent  freshwater  heads  as 
represented  in  piezometer  1,  and  therefore  the  environmental  heads  which  actually 
represent  the  level  to  which  the  saltwater  will  rise  as  shown  in  piezometer  2 have  to  be 
calculated.  These  values  will  be  used  to  interpret  results  and  calibrate  models. 
Eliminating  the  pressure  terras  and  soivtog  for  heads,  the  conversion  equations  arc 
derived  from  equations  (2-2)  and  (2-3)  to  be: 


Considering  a representative  elementary  volume  in  a porous  medium  (Figure  2-4). 
the  mathcraalical  expression  fijrthe  conservation  of  mass  is  written  os: 

( 3(/3p)'l  djpS) 

I & ^ & J ■ ft  ' 

whereyiis  the  fluid  density  (ML'’],  q is  the  specific  discharge  vector  [LT'],  p is  the 
density  of  water  entering  front  a source  or  leaving  through  a sink  [ML'’],  q,  is  the 
volumetric  flow  rate  per  unit  volume  of  aquifer  representing  sources  and  sinks  [T'],  ^is 
porosity  [dimensionless],  and  r is  time  [T]. 
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-if 

-if 


<2-41) 


and  finally,  Ihe  governing  equoiion  for  variable-density  (low  in  lernis  offrcshwalerhcad 
as  used  in  SEA  WAT  becomes  (Guo  and  Langevin.  2002): 


(2^3) 


where  the  subscripts  aand  ^represent  the  principal  directions  of  permeability  parallel  to 
the  bedding,  yrepresenu  the  direction  normal  to  the  bedding,  "nie  terras  ka,  kj  and  k,  are 
the  permeabilities  in  these  directions  (L’j,  and  5,,  Sp,  and  5,  are  the  angles  between  the 
respective  coordinate  axes  and  the  upward  vertical  direction,  as  described  in  the 
documentation  of  SEA  WAT. 

Variabie  density  groundwater  flow  causes  the  redistribution  of  solute 
concentration,  which  in  turn  alters  the  density  field.  This  redistribution  significantly 
affects  the  groundwater  movement.  Therefore,  it  is  important  to  describe  the  solute 
transport  process,  couple  it  with  the  flow  equation,  and  solve  the  resulting  equation. 

The  transport  of  solute  mass  in  groundwater  flow  is  brought  about  by  adveclioti, 
diffusion,  and  dispersion.  The  three-dimensional  equation  that  describes  this  process  is 
described  as  (Zheng  and  Bennett.  I99S): 


(k-i,  ...J9)is  the  rate  of  solute  production  or  decay  in  reaction  k of  N different  reactions 


^ = V-(DVO-V-{vC)-^C.*'£r, 


(2-44) 


where  D is  the  h 


nt  [L^'l,  P is  the  fluid  velocity  (LT 


],  Ct  is  the  solute  concentration  of  water  entering  from  sources  or  sinks  (ML*^],  and  Ri, 


SEA  WAT  solves  eoncemralion  in  leims  of  loial  dissolved  solids  (TDS).  The 
fonxiulation  of  a density  lerni  Ifoui  the  conceolsation  is  based  on  an  empirical  relationship 
(Baxter  and  Wallace.  1916): 

P = Pr<-SC  (2^5) 

where  the  term  £ represcnls  a dimensionless  constant  with  on  approximate  value  of  0.7 
for  salt  concentrations  ranging  from  aero  to  that  of  seawater  (roughly  3S.OOO  ppm). 

The  term  C represents  the  salt  concentration  |ML'’|.  When  the  empirical  equation 
of  Baxter  and  Wallace  (1916)  is  differentiated  with  respect  to  C.  the  resulting  relationship 
becomes: 


ec  ' 

When  the  term  that  relates  density  and  concentration  is  substituted  into  the  variable 
density  flowequadon  of  SEA  WAT,  the  governing  flow  equation  becomes: 


da 


1)+— 


-p,  az, 


dh,  Of 


The  linear  formulation  of  Baxter  and  Wallace  (1916)  is  considered  applicable  for 
normal  seawater  concentration  ranges.  However,  if  the  concentration  is  much  higher  or 
the  composition  of  seawater  is  different,  then  another  more  appropriate  formulation 
would  have  to  be  developed.  The  governing  partial  differential  equation  is  numerically 
solved  using  die  appropriate  boundary  and  ioitial  conditions. 


CHAPTER  3 

SEARCH  AND  OPTIMIZATION  IN  WATER  RESOURCES  MANAGEMENT 
OplimiziitJoa  Coonpu  and  Methods 
Opliniisstjon  modeling  is  pan  of  a general  class  of  problems  that  require 
minimization  or  meximization  of  an  objective  function  subject  to  a set  of  constraints. 
There  arc  several  types  of  optimization  modeling  based  on  the  objective  function  and  the 
constraints  to  be  imposed  (Walker  et  al.,  I99S).  The  mathematical  techniques  of  search 
and  optimization  provide  a formal  means  to  make  decisions  on  problems  that  require 
considering  various  options  or  best  alternatives  (Spall,  2003).  In  general,  search  and 
optimization  methods  are  used  to  obtain  iocal  or  global  optimal  values  or  a set  of  feasible 
solutions  for  a given  problem.  Global  optimization  is  a search  process  where  the 
objective  is  Rnding  the  best  possible  set  of  paramelcis  to  optimize  an  objective  function. 
Global  optimization  problems  generally  fall  within  the  broader  class  of  nonlinear 
programming.  There  are  two  groups  of  optimization  pmbicms.  combinoiorial  (or 
discrete)  and  continuous  variable  problems  (Colleilcand  Siatiy,  2003).  Search  methods 
currently  in  use  can  broadly  be  classified  as  caJculus-baacd.  cnumcrative,  and  random 
(Goldberg,  1989). 

Cslculus.based  search  methods  search  for  the  optimal  value  of  an  objective 
ftmclion  in  two  ways,  either  directly  or  indirectly  (Goldberg,  1989).  The  indirect  way 
requires  finding  local  extrema  by  solving  a set  of  equations  that  is  usually  non-linear  by 
selling  the  gradient  of  the  objective  (unction  to  zero,  The  search  in  this  case  is  restricted 
to  points  where  the  slope  is  zero  in  all  ditoctions.  The  direct  approach,  on  the  other  hand, 
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searches  for  local  oplima  by  picking  points  on  ihc  fiinclionand  moving  in  a direction 
relaled  to  the  local  gradienl-  In  other  words,  this  is  a hill  climbing  lechnique  where  Ihe 
local  heal  is  obtained  by  climbing  the  funclion  in  Ihc  steepest  permissible  direction. 
Because  of  their  local  scope  in  Ihe  gradienl  analysis,  calculus-based  techniques  yield 
local  optima.  Moreover,  oblaining  gradients  and  derivatives  is  not  always  possible.  For 
these  reasons,  the  general  class  ofcalculus-bascd  methods  lacks  robustness. 

Enumeralive  methods  work  within  a finite  search  space  or  a discretized  infinite 
search  space  to  search  for  the  optimal  value  (Goldberg,  I9S9).  The  algorithm  considers 
the  objective  funclion  values  at  every  point  in  ^acc  one  at  a time  One  such  example  is 
dynamic  programming.  Due  lo  the  requirement  of  function  value  evaluation  at  every 
point  in  space,  (wholher  ilis  lirule  search  space  or  discretized  infinite  space)  enumeralive 
techniques  arc  not  very  efiicienL 

Random  techniques  try  to  obtain  the  global  optima  by  randomly  searching  and 
saving  a SCI  of  optimal  values  and  coniinuing  the  somewhat  slochaslic  process  using  an 
algorithm  uniil  the  slopping  criterion  is  met  (Goldberg,  1989).  The  search  docs  not 
require  dirfcremialing  the  objective  Function  to  obtain  gradients,  thus  increasing  the 
suilablliiy  and  applicability  of  Random  techniques.  Overall,  Random  search  techniques 
are  computationally  intensive  and  fall  short  with  regard  lo  elficiencydue  to  Ihe  inherent 
naniieofthc  search  method.  With  current  high  speed  and  large  memory  coropulers  the 
efficiency  problem  could  be  minimized. 

Through  the  years  numerous  techniques  of  optimization  have  been  developed  lo 
solve  a general  class  of  optimization  problems.  Some  of  Ihc  major  techniques  of 
optimization  are  described  very  briefly  in  Ihe  following  paragraphs. 


5! 

Discrele  Opiimizolion 

Discrete  opiimizalion  Is  used  to  describe  a category  of  optimization  methods  where 
the  best  alleniaUve  is  selected  from  a finite  set  of  feasible  solutions  as  derined  in  the 
objective  function.  The  variables  used  in  the  objective  function  of  discrete  optimization 
are  restricled  to  assume  only  integer  values.  The  most  common  techniques  used  in 
solving  discrete  optimization  problems  included  integer  programming,  heuristic  search 
techniques  such  as  genetic  algorithms,  simulated  annealing,  and  artificial  neural 
networks,  and  global  search  techniques 
Continuous  Optimization 

Continuous  optimization  refers  to  a general  class  of  optimization  problems  in 
which  all  Ihe  variables  ace  allowed  to  take  values  from  subintcrvals  of  the  real  line.  Thus, 
continuous  opiimization  allows  decision  variables  delined  in  the  objective  function  to 
have  real  orcontimious  values. 

The  major  optimization  methods  used  in  both  discrete  and  continuous  optimization 
are  briefly  described  in  the  following  paragraphs. 

Linear  Programming.  Linear  programming  is  based  on  linear  objective  luitclion 
and  constraints  with  continuous  decision  variables.  Theolqective  ofan  LP  problem  Is  to 
find  the  minimum  or  maximum  of  a linear  objective  function  subject  to  linear  constraints. 

Integer  Programming.  Integer  programming  Is  a technique  that  utilizes  linear 
objective  function  and  constraints  with  integer  decision  variables. 

Mixed  Integer  Programming.  The  optimization  technique  of  mixed  integer 
programming  was  developed  to  solve  linear  objective  functions  with  linear  constraints 
andinteger  and  continuous  decision  variables.  MIP  Is  an  extension  of  integer 
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prognunmuig  where  Ihe  problem  has  one  or  more  Integer,  real,  or  binary  ilecision 
vahables 

Dynamic  Programming.  Dynamic  programming  iso  technique  used  to  optimize 
mullislage  processes.  The  applicability  of  dynamic  programming  in  water  resources 
optimization  problems  was  promoted  by  Ihe  capability  of  dynamic  programming  to 
translate  nonlinear  and  slochastic  problems  into  adynamic  programming  formulation. 
This  was  achieved  by  decomposing  complot  problems  with  large  number  of  variables 
into  a series  of  sub  problems  which  are  solved  ileralively. 

The  primary  advantage  of  dynamic  programming  technique  ties  in  its  capability  to 
handle  discrete  variables  and  nonconvex.  noncontinuous,  and  nondiffercniisblc  Functions. 
It  can  also  lake  into  account  alochasilc  variability  by  a simple  modificalion  of  Ihe 
dclciminisijc  procedure.  In  order  to  obtain  a solution,  dynamic  programming  requires  the 
objective  function  to  satisfy  conditions  of  separability  and  monolonicity.  An  objective 
function  is  said  to  be  sepatubie  if  decomposed  components  of  the  objective  function  are 
independent  of  each  other.  Similarly,  an  objective  function  is  said  to  be  monoionic  if 
improvements  in  the  decomposed  individual  components  of  the  objective  function  lead  to 
improvements  in  the  objective  function.  The  major  drawback  of  dynamic  programming 
Is  dimensionality  due  to  the  large  number  of  possible  solution  designs.  The  size  is 
generally  an  exponential  function  of  locations  and  periods.  However,  differential 
dynamic  programming  overcomes  the  problem  of  dimensionality  by  eliminating  the  need 
for  discretization  of  Ihe  control  and  slate  vector,  and  by  introducing  stepwise 
decomposition.  Difrcrenlial  dynamic  programming  is  characterized  by  a linear  growth  in 
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cocnpullng  effon  with  resp«cl  lo  the  mimber  of  stages  or  planning  periods,  and  it  also 
exhibils  quadratic  convergence. 

Nonlinear  ProgrammlDg.  When  both  the  objective  function  and  the  decision 
variables  are  nonlinear,  the  problem  can  be  salved  using  the  techniques  of  nonlinear 
programming.  Many  groundwater  planning  and  management  models  involve  nonlinear 
objective  Aincllons  and  eonsiroinis.  The  nonlinearilies  could  be  due  lo  nonlinear  coal 
functions,  nonlinear  equations  governing  the  flow,  panicularly  for  unconfined  aquifers, 
nonlinearities  in  the  governing  equations  for  solute  transport  in  groundwater,  other  types 
ofnonllnear  physical,  and  managerial  objective  functions  and  constraints. 

Combinatorial  optimization  deals  with  problems  that  have  a linear  or  nonlinear 
function  defined  over  a finite  but  large  set  ofproblem  space.  Examples  of  combinatorial 
optimization  problems  include  network  problems,  scheduling,  and  transportation,  hurilily 
location  and  layout,  and  optimal  design  of  waterways  and  bridges.  If  the  function  is 
piecewise  linear,  the  combinatorial  problem  can  be  solved  exactly  with  a mixed  integer 
program  method,  which  uses  branch  and  bound.  Combinuiorial  optimization  problems 
generally  consist  of  finding  the  optima  for  a function  with  constrainu  and  encoded  as 
integer  vectors  or  binary  variables.  The  methods  used  to  solve  combinatorisl 
optimization  problems  include  branch  and  bound,  fixed  heuristics,  and  mcla  heuristics. 

In  the  general  class  of  heuristic  methods,  evolutionary  algorithms  such  as  simulated 
annealing,  tabu  search,  and  genetic  algorithms  have  been  successfully  used  to  solve 
combinatorial  optimization  pioblems. 

Evolutionary  Algorithms 

Evolutionaiy  algorilhms  ore  a new  class  of  heuristic  search  and  optimization 
methods.  These  melhods  have  been  called  global  optimization  methods  or  gradient  free 
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methods.  The  term  evolutionary  algoriUtms  has  also  been  used  to  describe  some  of  them 
since  they  are  derived  from  natural,  biological  principles.  Evolutionary  metlrods  and 
artificial  intelligence  have  been  gaining  increasing  attention  in  the  past  few  years. 
Evolutionary  mclhods  try  to  emulate  the  evolution  of  life  through  natural  selection  artd 
use  similar  procedures  to  obtain  the  best  fining  optimal  solution  for  a given  problem. 
Anificlal  inlelligenee  methods  on  the  other  hand  try  to  train  the  system  through 
exonrples,  in  a fashion  similar  to  the  way  a child  learns  about  the  environment.  Then, 
once  a sumdenl  knowledge  database  has  been  acquired  through  training,  the  system  is 
required  to  recognize  the  best  solution  or  simulate  outputs  for  a given  set  of  inputs. 
Genetic  algorithms  (GAs)  are  one  of  the  popular  evolutionary  algorithms  used.  Artificial 
neural  networks  (ANNs)  on  the  other  hand  have  become  increasingly  applicable  as 
artificial  intelligence  methods  to  perfonn  simulation  and  optimization  in  a number  of 
fields. 

Simulnted  anncaliDg 

Simulated  annealing  is  a process  based  on  the  principles  of  thermodynamics,  the 
basis  of  which  is  a process  in  which  the  temperature  of  a solid  is  raised  lo  bring  its 
molecules  into  a highly  mobile  stale  and  then  cooled  lo  form  a Iow*energy  crystalline 
structure.  Simulated  annealing  was  introduced  by  Metropolis  et  ah  (1953)  and  is  used  to 
approximate  the  solution  of  very  large  combinatorial  optimization  problems  (Kirkpatrick, 
1983).  The  technique  origtnales  from  the  theory  of  slalistioal  mechanics  and  is  based 
upon  the  analogy  between  the  annealing  of  solids  and  solving  optimization  problems. 

The  four  key  ingredients  for  the  ImplemenUlion  ofslmulaled  annealing  ate:  the 
definition  of  configurations,  a generation  mechanism,  i.e.,  the  definition  of  a 
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ncighbcrhood  on  the  coniigiirolion  spoce,  the  choice  of  a cosi-funclion,  and  a “cooUng" 
schedule. 

Tabu  search 

Tabu  search  is  a combinalorial  oplimization  mcihod,  analo|ous  lo  gradieni-dcsccnt 
search  vsiih  memory,  and  ii  is  designed  to  avoid  local  minima  (Clover  1986, 1990).  The 
memory  preserves  a number  of  previously  visited  slates  along  with  a number  of  stales 
that  might  be  considered  unwanted.  This  information  is  stored  in  a tabu  list.  The 
definition  of  a state,  the  area  around  It  and  the  length  of  the  tabu  list  are  critical  design 
parameters,  in  addition  to  these  tabu  parameters,  two  extra  parameters  are  often  used: 
aspiration  and  diversificaiian.  Aspiration  is  used  when  all  the  neighboring  states  of  the 
current  state  are  also  included  in  the  tabu  list.  In  that  case,  the  Ubu  obstacle  is  ovemdden 
by  selecting  a new  stale.  Diveisificalion  adds  randomness  to  this  otherwise  detemtlnisHc 
search,  if  the  tabu  search  is  not  converging,  the  search  Is  reset  randomly. 

Genetic  algorithms 

Genetic  algorithms  are  search  algorithms  that  are  based  on  the  principles  of 
evolution.  Therefore,  the  formulation  of  a genetic  algorithm  involves  natural  selection, 
inheritance,  and  mutation.  Through  these  processes,  OAs  utilize  survival  of  the  fittest 
among  string'  structures  with  a structured  yet  randomized  information  exchange  to  form 
a search  algorithm  (Goldberg,  1989).  The  development  of  GAs  as  search  and 
optimization  tools  was  fust  pioneered  by  Holland,  at  the  University  ofMIchigan  In  1975 
(Goldberg.  1989).  The  driving  force  behind  the  application  of  GAs  was  robusmess. 
which  is  the  balance  between  efficiency  and  efficacy  (Goldberg,  1939). 


'Suuigs  are  a set  ofonifieial  cieaiures  rormed  tram  bus  and  pieces  of  the  finest  or  The  old  Bencrsnon. 
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Genclic  algoriUims  dilTer  from  ihc  Iraditional  oplimizaiion  methods  in  many  ways. 
The  following  are  unique  features  ofGAs  (Tang  and  Mays,  1998.  and  Goldberg,  1989); 

• Ability  to  use  the  objective  itself;  not  ils  derivatives; 

• Ability  to  search  from  a population  of  decision  variable  sets,  not  a single  variable 

• Inherent  random  property  that  helps  avoid  local  optima;  and 

• The  use  of  probabilistic  transition  rules  as  compared  to  deterministic  rules. 

As  a result,  GAs  are  very  suited  to  handle  convex  and  non-linear  problems  by 

ensuring  a global  or  near  global  optimal  solution. 

There  are  five  components  in  the  development  of  a GA  solution  for  optimization 
problems  (Adeli  and  Hung.  1995).  These  steps  are  encoding,  initialization  of  the 
population,  fitness  evaluation,  evolution  performance,  and  working  parameters: 
Encoding:  Encoding  Involves  generating  a set  of  decision  variables  of  the  problem  and 
encoding  it  into  a string  called  a chromosome.  The  most  common  method  uses  a binary 
siring  to  represent  this  vector  of  decision  variables.  Sometimes,  integer  and  real  coding 
are  also  utilized  instead  of  binary  coding.  If  there  ore  m decision  vanables  in  an 
optimization  problem  and  each  variable  is  encoded  as  a binary  string  of  length  a,  Ihcn  the 
resulting  chromosome  will  be  an  n x m string. 

Populaiian  Imiialaalioa:  The  initiabzation  is  generally  performed  by  random 
number  generation  by  assigning  an  equal  probability  of  selection  to  each  bit. 

Fitness  Evalualion-.'ne  fitness  evaluation  helps  determine  the  probability  that  a 
chromosome  will  be  selected  as  a parent  chromosome  to  generate  new  chromosomes. 
This  is  achieved  through  the  use  of  a scaling  function  8.  which  maps  objective  fiinction 
values  to  an  appropriate  positive  fitness  value. 
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Evolulion  Performance-.  Generally  Ihere  are  ihree  operalors  lo  peribnn  evolution 
between  two  consecutive  chromosome  populations.  These  are  selection,  crossover  and 

Working  Parameun.  The  development  of  working  parameters  requires  predefining 
a set  of  parameters  to  guide  the  GA.  The  parameters  are  chromosome  length,  population 
siae,  crossover  rale,  mutation  rate,  and  stopping  criterion. 

Application  of  Optimization  in  Groundwater  Management 

There  are  a number  ofmaiheraalical  tools  used  lo  optimize  the  utilization  of  water 
resources  subject  to  various  constraints.  The  constraints  in  most  cases  require  meeting 
some  ecological  needs,  quantitative  and  qualitative  hydrologic  constraints,  forecasted 
demands,  and  in  some  cases  risks. 

The  formulation  of  practical  water  resources  optimization  problems  is  characterized 
by  convexity  and  non-linearity.  As  a result,  obtaining  second-order,  and  even  In  some 
cases,  firsl-order,  derivatives  of  the  objective  fiinclions  that  are  required  by  many 
iradilionaJ  optimization  methods  is  dilficuli  if  not  impossible  (Tang  and  Mays.  I99S). 
Moreover,  it  ia  difficult  to  ensure  lhai  the  solution  oblalned  is  the  globally  optimal 
solution  for  a given  problem.  A number  of  optimization  models  have  been  developed  to 
address  the  issue  of  groundwater  manogcmcni.  Most  models  are  developed  to  foimulalc 
optimal  pumping  strategies  by  employing  cost  minimization  or  pumping  maxiiuizaiion 
subject  lo  water  quality  and  quantity  conslniinis. 

Oplimizalion  Coupling  Methods  In  Groundwater  Systems 

Optimization  algorithms  developed  for  management  of  groundwater  systems  have 
to  be  coupled  with  a simulation  model.  The  coupling  of  single-objeclive.  linear 
groundwater  optimization  models  is  broadly  classified  on  the  basis  oflhe  approaches 


used.  The  two  approaches  ate  embedded  and  response  matrix  approaches  (Gorelick, 
1983). 

Embedding  methods  involve  development  of  a linear  set  of  equations  known  as 
response  equations  from  the  partial  differenllaJ  equalions  of  groundwater  How  and 
transport.  The  solution  of  the  resulting  equalions  is  obtained  by  using  numerical  methods 
such  as  riniiedifrerencc  or  finite  elements.  Applications  of  this  approach  have  been 
successfully  tested  by  Aguado  and  Remson  (1974)  and  Willis  and  Liu  (1984). 

The  response  matrix  of  groundwater  optimization  problems  is  developed  from 
hydraulic  orwalet-quality  models  of  the  aquifer  system.  The  simulation  model  is  used  to 
generale  the  response  oflhe  aquifer  system  to  pumping  or  recharge.  The  resulting  matrix 
ofcoemcienis  is  used  to  convert  pumping  stresses  into  linearized  changes  In  drawdowns 
and/or  changes  in  concentration  of  the  ambient  water.  The  set  of  unit  responses  is 
formulated  for  individual  wells  and  superimposed  to  model  multiple  wells,  which  form  a 
set  of  linear  equalions.  The  resulting  set  of  equations  is  solved  systematically.  However, 
if  the  problem  is  formulated  as  a non-linear  problem,  then  the  solution  is  obtained 
iteralively.  Examples  of  response  matrix  application  in  multi-objeclivc  optimization 
include  (Louie  el.  ah.  1 984)  and  (Yeh  eu  al.,  1992). 

Linear  programming  techniques  have  been  used  in  groundwater  opiimizalion 
problems  where  both  the  objective  function  and  constraints  were  linear  flinctions- 
Aguado  eL  al.  (1974)  applied  linear  programming  to  predict  optimum  number  and 
locations  of  wells,  and  rales  of  pumping  needed  to  maintain  ground  water  levels  below 
specified  elevations  undm  steady  state  conditions.  The  objective  was  to  minimize  total 
pumping  while  mainlaining  steady  state  groundwater  hcods  below  some  assigned  value 


during  dewatering  of  a large  dry  dock  excavuiion  area.  Goielick  and  Remeon  (19S2) 
incorporated  the  steady  state  finite  diiference  form  of  the  solute  transport  equation  as  an 
embedded  constraint  to  maximize  waste  disposal  at  two  locations  while  protecting  water 
quality  at  supply  wells  and  maintaining  an  existing  waste  disposal  facility.  The  technique 
has  been  used  extensively  in  groundwater  optimization  pioblems.  Linear  programming 
based  management  models  were  developed  by  embedding  the  governing  di/TerenliaJ 
equations  as  constraints  in  the  formulation  (Aguado  and  Remson,  1974),  The  objective 
hmetionwas  formulated  as  maximization  oflhe  hydraulic  heads  at  specific  locations 
subjecl  10  sum  ofpioduciioa  mie  and  non  decreasing  head  directional  monolonicity 


Mixed-integer  programming  has  been  used  to  solve  optimization  problems  with 
linear  objective  function  and  linear  constiainls  in  which  some  of  the  variables  can  lake 
only  integer  valuea.  These  types  of  requirements  arise  in  situalions  where  the  decision 
variables  oflhe  management  problem  require  a logical  decision  prompt  such  as  yes  onto, 
or  in  situalions  where  the  decision  variables  decide  parameters  such  as  the  number  of 
installations  and  locations,  in  instances  where  all  the  variables  lake  only  integer  values, 
the  problem  invariably  reduces  lo  on  mieger  programming  problem.  A mixed-integer 
programming  with  the  response  matrix  approach  was  used  to  determine  an  optimal 
dewatering  scheme  in  the  foundation  design  of  an  eicctronuclear  piani  (Galeati  and 
Gambolati,  1988).  The  response  matrix  was  generoled  using  a three-dimensional  finite 
element  model  for  steady  state  flow  conditions.  The  best  locations  for  a specified 
number  of  pumping  wells  were  determined  by  using  the  branch  and  bound  method  to 
solve  the  MIP  objective  function  (Rosenwald  and  Green,  1974).  A planning  model  for 
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ihe  oplimol  conjunctive  use  of  groundwatet  and  surface  water  resources  was  formulated 
and  applied  using  MIP  (Willie.  1976).  Olher  applications  of  MIP  in  groundwater 
management  are  presented  in  Rosertwald  and  Green  (1974)  and  Aguado  and  Remson 
(1980).  Similar  applications  include  (Aguado  and  Remson,  1980).  (Misirli  and  Yazicigil. 
1997),  and  (Datto  and  Dhiman,  1995). 

Most  groundwater  optimization  problems,  especially  imeonfined  aijuifers.  exhibil 
nonlinearity.  There  arc  a number  of  factors  that  make  the  groundwater  system  nonlinear. 
These  factors  include  nonlinear  cost  funciions,  nonlinear  flow  and  Iranspon  equations, 
and  nonlinear  objective  functions  and  constraims.  Nonlinear  programming  has  been  used 
lo  solve  nonlinear  opiimization  problems  in  water  resources  management.  Gordick  el.al. 
(1984)  developed  a planning  model  lo  determine  the  oplimal  design  of  reclamation 
schemes  for  contaminated  groundwater  systems.  The  model  combined  a nonlinear, 
distributed  parameter  groundwater  flow  and  solute  transport  simulation  model  (SUTRA) 
with  a nonlinear  optimization  model  (MINOS).  Colarulloel.  al.  (1984)  developed  a 
hydraulic  management  and  plume  stabilization  model  for  a partially  contaminated  aquifer 
used  for  both  water  supply  and  waste  disposal.  Their  model  consisted  of  a quadratic 
objective  (unction  with  groundwaler  velocity  constraints. 

The  exponential  form  of  fixed  charges  (installation  cost)  was  incorporated  into  Ihe 
objeclive  ftmetion  end  the  problem  was  solved  as  a concave  minimization  problem 
(Karatzas  and  Finder,  1993).  They  applied  the  outer  approximation  method  to  concave 
global  minimization  problem  over  a convex  corapacl  set  of  constrainu.  The  cost  analysis 
of  an  aquifer  restoraiion  system  was  also  performed  using  quadratic  programming 
(Lelkoffand  Gorelick,  1986).  The  design  and  cost  analysis  of  aquifer  restoraiion  systems 


were  formulaled  by  sssmning  adveclive  groiindwaier  How  and  ineorporaling  mardaiion 
of  solutes  due  losorplion.  Olher  applications  of  nonlinear  oplimiaation  include  Wang 
and  Ahifeld  (1994),  and  Hallaji  and  Yazicigil  (19%). 

The  objective  function  defined  in  terms  of  minimizing  the  present  value  of 
pumping  costs  was  solved  using  quadratic  programming  combined  with  an  algebraic 
technological  function  (Maddock,  1972a).  These  technological  functions  or  response 
coefficients  were  defined  as  the  changes  in  drawdown  induced  by  unit  pumping  at  each 
well.  Maddock  (1972b)  also  used  quadratic  programming  with  a nonlinear  objective 
function  and  linear  constraints  to  solve  the  problem  of  management  of  an  unconfined 
aquifer.  Mixed-integer  quadratic  programming  was  used  lo  minimize  the  pumping  costs 
plus  fixed  cosis  for  well  and  pipeline  coostructlon. 

Dynamic  programming  provides  tools  to  solve  multistage  decision  problems  (Das 
and  Dana,  2001).  Solving  multistage  problems  using  dynamic  programming  techniques 
requites  decomposing  the  multistage  problems  into  single  stage  problems.  Individual 
single-stage  problems  may  then  be  solved  by  any  method  ofqitimization.  The  advantage 
of  using  the  dynamic  programming  is  its  capability  lo  deal  with  non-convex,  non- 
cominuous  and  non-diffetentiable  fiinclions.  and  discrete  variables,  A penally  function 
approach  was  used  in  conjunction  with  a dynamic  programming  technique  lo  obtain 
solutions  to  the  constrained  optimal  eonnol  problem  with  a large  number  of  constraints 

(Chang  el  al..l  992).  A hyperbolic  penalty  ftinction  was  used  to  incorporate  the 
conslcainls.  A lime  varying  optimization  of  groundwater  remediation  with  variable 
management  and  simulaiion  periods  was  also  solved  using  dinerential  dynamic 
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programming  (Culver  and  Shoemaker,  1992).  Yakowitz  (1982)  has  given  a detailed 
roview  of  various  DP  approaches  used  in  water  resouiccs  management. 

A difTerenliai  dynamic  programming  algorithm  was  used  to  solve  unsteady 
nonlinear  groundwater  management  problem  where  the  management  problem  was 
formulated  as  an  optimal  control  problem  (Jones  CL  aJ.,  1987).  Culver  and  Shoemaker 
(1992)  presented  a differential  dynamic  programming  algorithm  for  time  varying 
optimization  of  groundwater  remediation  in  which  management  periods  were  different 
from  the  simulation  periods. 

Combinatorial  optimization  techniques  search  for  optimal  values  within  a well 
defined  discrete  problem  space.  Evolutionaiy  methods  generally  fall  under  the  heuristic 
methods.  Application  of  evolutionary  methods  in  groundwater  optimization  problems 
have  been  researched  and  documented  in  the  past  two  decades.  Some  examples  include 
Zheng  and  Wang  (1996),  Johnson  and  Rogers  (2000),  Mairyott  eL  al-(1993).  Morshed 
and  Kaluarachchi  (1998),  and  Cheng  el  lU.  (2000). 

Formulation  of  Multi-objective  Optimization  Problems 

Groundwater  management  problems  ollen  require  considering  moro  than  one.  and 
usually  conflicting,  managemeol  objectives.  For  example,  a management  problem  where 
it  is  required  to  maximize  the  amount  of  water  to  be  pumped  from  an  aquifer  and 
minimize  the  drawdown  represents  a problem  with  confUeting  management  objectives. 
Due  to  the  nature  of  the  management  issues  involved,  groundwaicr  management 
problems  are  often  formulated  as  multi-objective  optimization  problems.  Multi-objective 
groundwater  management  models  are  often  formulated  as  mathematical  programming 


problems  wiih  mony  conflicting  objectives  that  require  u set  of  Pareto  optimal . non- 
dominaled,  or  non-inferior  (efTicicnl)  solutions  (Das  and  Datta,  2001). 

The  mathematical  definition  of  global  optimality  or  Pareto  optimality  for  a vector 
of  decision  variables  i'  defined  in  a set  of  solutions  C is  given  as  follows; 

A vector  s C is  said  to  be  (globally)  Pareto  optimal,  or  a (globally)  efficient 
solution,  or  a non-dominated,  or  a non-inferior  point  for  multi-objective  opUmization 
problem  (MOP)  ifand  only  ifthcrc  is  no  js  C such  that  the  functional  value  of  the 
objective  function /,(.t)S/(a’)for  all  is  (1.2,. ,.n(  with  at  least  one  strict  inequality’. 
Multi-cnteria  optimization  methods  typically  assume  that  a multi-criteria  problem  is 
converted  into  an  auxiliary  parametric  single-objective  problem  whose  solution  provides 
a Pareto-optimal  point.  The  Pareto  optimal  solutions  are  represented  by  a trade-off  curve 
relating  two  conflicting  objectives  (Das  and  Datta.  2001).  Multiple  objective 
management  models  generally  attempt  to  develop  such  trade-off  curves.  Multi-objective 
optimization  problems  arc  mostly  solved  by  combining  the  multiple  objectives  into  one 
scalar  objective  whose  solution  is  a Pareto  optimal  point  for  the  original  multi-objective 
problem. 

A multi-objective  groundwater  quality  manogemenl  model  was  developed  for  a 
regional  unconfined  aquifer  (Willis,  1977).  The  multi-objective  model  was  formulated  to 
minimize  the  operational  cost  of  pumping,  minimize  the  minimum  hydraulic  head  within 
the  injection  zone  of  the  aquifer,  and  maximize  the  injection  rale.  The  response  equations 

' Asolulion  is  called  PareloAiounBKof  elficieat)  solctiaii,  iflbereisDOolhersolurionforwhichat  least 
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were  used  as  embedded  conslraints  in  the  multi-objeclive  oplimlzalion  model.  The 
optimization  model  was  bounded  by  constralms  of  water  demand  and  waste  load,  storage 
capacity,  maximum  pumping,  maximum  and  minimum  permissible  heads,  response 
equations,  and  non-negativity  of  the  state  and  decision  variables. 

A multi-objective  linear  programming,  coaslal  aquifer  management  model  was 
developed  to  manage  the  problem  of  saltwater  intrusion  (Shamir  el  aJ.,  1984).  The 
management  model  had  four  objective  functions  which  included,  groundwater  contour, 
location  of  the  toe  of  the  seawater-freshwater  interface,  concentration  distribution,  and 
minimum  energy  for  pumping  and  recharge.  An  approximate  linearized  expression  was 
used  to  incorporate  the  location  of  the  seawater-freshwater  interface,  and  a constraint 
method  ormulii-objeotive  analysis  (Cohon  and  Marks  1978)  was  used  to  obtain  trade-off 
functions  between  paiia  of  objectives. 

Muld-objective  programming  methods  were  used  to  determine  the  optimal 
groundwater  management  schemes  in  a hypothetical  multi-aquifer  system  (Yazicigil  and 
Rashccduddin,  1987).  The  management  problem  wos  solved  by  using  a combination  of 
embedding  technique  and  linear  programming.  A single  objective  case  was  applied  to  a 
transient  aquifer  management  while  the  multi-objective  model  was  developed  fora 
steady  stale  aquifer  management.  Trade-off  curves  were  developed  for  hydraulic  heads 
and  water  withdrawal  targets  using  consuainis  and  weighting  methods.  Similarly,  a 
multi-objective  management  model  was  used  to  determine  the  optimal  management  of  a 
muln-layer  aquifer  system  (Yazicigil,  1990). 


CHAPTER  4 

DEVELOPMENT  OF  A GEOGRAPHIC  INFORMATION  SYSTEMS  MODEL  TO 
INVESTIGATE  THE  IMPACTS  OF  SEA  LEVEL  RISE 

Concepts  ond  Methods 

Gcogcaphiclnfoiinalion  Systems  (CIS)  provide  the  tools  and  methods  required  to 
perform  spmiiil  analysis  on  georeferenced  data.  Using  GIS  feciliialcs  investigation  of  the 
impacts  of  sea  level  rise.  The  visual  dimension  provided  by  GIS  makes  it  possibie  to 
better  understand  the  changes  in  sea  level  with  respect  to  space  and  time  for  a given  area. 
The  primary  interest  in  such  analyses  is  to  determine  how  far  seashore  would  move 
inland  if  the  sea  level  were  to  rise  by  a cenain  magnitude.  Subsequent  investigations 
could  focus  on  Identi5'ing  and  delineating  areas  that  would  be  affected  as  a result  of  sea 
level  rise,  assessing  the  landuse  and  landcover  that  would  be  impacted,  and  quantifying 
the  loss  of  habiut  due  to  submergence.  Application  of  GIS  to  model  components  of  the 
ecosystem  have  been  described  in  Singh  and  Fiorentino  (1996).  Tooley  and  Jelgersma 
(1992),  Baltaglin  et  al.  (1993),  Ooodchild  ct  al.  (1995),  and  Michener  et  al.  (1994). 

The  analysis  of  sea  level  rise  primarily  relies  on  elevation  data  and  digital  elevation 
data  processing  techniques.  Digital  Elevation  Models  (DEMs)  have  been  used 
extensively  to  analyze  terrains  and  derive  design  parameters.  Elevation  data  in  DEMs  is 
stored  in  raster  fdrmaL  which  is  one  of  the  data  formats  currently  in  use.  Raster  GIS  deals 
with  data  stored  and  manipulated  on  a cell  by  cell  basis.  The  DEM  of  an  area  therefore, 
consists  of  cells  where  the  required  data  are  slored  in  a reco^zable  format  The 
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developers  and  solution  providers.  AreView  is  one  of  the  desktop  CIS  aitalysis  softwares 
developed  by  ESRl.  An  AreView  extension  called  Spatial  Analyst  make*  it  possible  to 
perform  spatial  analysis  on  Roster  data. 

In  order  10  create  maps  ofareas  vulnerable  to  sea  level  rise  using  seleeled 
scenarios,  the  simulation  requires  analyzing  bpul  values  on  a cell  by  cell  basis.  The  cell 
by  cell  analysis  is  a compulationally  intensive  process  that  AreView  is  not  well  equipped 
to  handle  as  a stand-alone  software.  Asa  result,  the  core  compulations  have  to  be 
performed  using  a separate  tool  that  Is  compulalionally  robusi  and  that  can  inlerface  with 
AreView. 

AreView  was  originally  wrinen  in  an  object  oriented  language  called  Avenue, 
although  currenlly  ESRl  is  shifting  to  the  more  widely  used  language  Visual  Basic.  The 
primary  lask  in  the  sea  level  rise  modeling  process  was  iherefore  customizing  AreView 
by  writing  source  codes  to  read  the  required  inpul  data  iheo  to  export  the  digital  elevation 
data  and  other  relevant  information  to  on  external  program  to  perform  the  computation, 
obtain  the  processed  data,  create  laslcr  files  of  the  result,  and  generate  maps. 

CIS  Model  Framework 

1he  CIS  model  developed  in  this  research  consists  of  on  AreView  inlerface  wrilten 
in  Avenue  and  a core  computational  algorithm  written  in  FORTRAN.  The  model 
currenlly  can  be  used  lo  investigate  impacts  of  sea  level  rise.  However,  any  spatially 
referenced  parameter  in  rasier  formal  could  be  studied  by  appropriate  modification  of  ihe 
corealgorilhm.  For  example  if  the  user  wants  lo  sludy  spatial  variability  of  rainfall,  Ihe 
core  algorithm  that  performs  the  analysis  is  modified  to  compute  Ihe  required  results 
while  the  main  inlerface  is  modified  to  read  the  required  raster  file  and  rainfall  data.  In 
Ihe  current  conllguralion.  the  input  data  to  Ihe  model  consists  of  the  DEM  of  an  area,  a 
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seiccicd  scenario  of  climaie  change,  a selecled  simul alien  period,  and  the  local  sea  level 
Irend.  The  IPCC  scenarios  provide  global  sea  level  rise  values  lhal  are  used  lo  calculalc 
ihc  [clallve  sea  level  rise  for  the  area.  The  ouipui  of  the  model  is  displayed  as  a map  lhal 
could  be  exported  lo  serve  as  a boundary  delinilion  map  for  a groundwater  flow  model. 
The  processes  involved  in  the  model  developed  as  pan  of  ihe  research  described  in  this 
dissertation  are  described  in  the  flow  chan  presented  in  Figure  4*1-  The  simulation 


model  is  named  CCGEN  (Climaie  Change  scenario  GENeralor), 


As  displayi^  in  Figure  4-1.  the  simulaiion  model  hoa  a modular  struclure  ihsl 
allows  further  development  and  addition  of  taster  based  analyas  packages-  CCGEN  was 
developed  and  tested  on  a peisonal  computer  running  Windows  2000  and  can  be  run  on 
any  computer  that  runs  Windows  versions  98  and  above,  a memory  (RAM)  of  at  least 
128  MB,  and  AreView  3x  with  Spatial  Analyst.  The  model  was  compiled  in  such  a way 
that  the  FORTRAN  code  is  embedded  within  the  Avenue  calling  program.  The  whole 
package  was  developed  as  an  Arc  View  project,  but  it  can  also  be  recompiled  as  an 
extension  for  ArcView3x  to  make  it  easily  portable. 

Algorithm  Development  and  Programming 

The  primary  factor  in  the  development  of  a sea  level  rise  modeling  algorithm  is 
chooslagasuiloblemalhefnalical  representation  that  defines  the  coastal  geomeuy  and 
processes.  Sea  level  rise  primarily  results  in  shoreline  retreat  due  to  the  combined  effects 
of  shore  erosion  and  inundation.  Inundation  is  the  result  of  continuous  and  progressive 
submergence  while  erosion  is  removal  of  materials  from  the  shores  and  beaches.  The  rise 
in  sea  level  and  the  subsequent  retreat  in  shoreline  are  presented  schematically  in  Figure 
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Figure  4-2,  Schematic  representation  of  shoreline  relreat  due  to  sea  level  rise. 


Changes  in  shoreline  profile  due  lo  sea  level  rise  could  be  quantified  using  a variety 
of  approaches.  The  simplest  method  in  use  is  based  on  the  concept  of  drowned  valley 
(Dean  el  a].,  1987).  in  which  the  preexisting  lopogtnphy  along  the  shore  is  fixed  and 
combined  with  sea  level  rise  lo  project  the  profile  of  the  new  shoteline  (Kana  et  ol., 

1 9S4).  Since  slope  is  the  cunlrolling  variable  in  this  approach,  coastal  areas  with  gentle 
slopes  will  be  subject  to  larger  flooding  and  shoreline  retreat  while  shores  with  sleep 
slopes  will  experience  little  horizomal  shoreline  displacement. 

There  are  a number  of  techniques  to  estimate  the  shoreline  retreat  primarily  due  to 
etosional  potential  of  sea  level  rise,  a few  of  the  methods  documented  in  (Dean  el  al., 

1 987)  are  described  as  fellows: 

1.  Historical  Trend  Analysis  (Lealheraian,  1984):  Exlrapolslion  of  historical  trends  or 


historical  trend  analysis 
method  accounts  for  varii 


calibiation  method  based  on  historical  tecords.  The 
ability  in  shoreline  response  to  coastal  processes, 
and  coastline  exposures. 


The  Bruun  Rule  (Bruun.  1962):  Bnnin's  rule  is  based  on  what  is  known  as  the 
equilibrium  beach  profile,  which  refers  to  a slalislical  average  profile  that  maintains 
its  form  apart  from  small  fluctuations,  including  seasonal  cfTects  al  a particular 
water  level  (Bruun,  1954).  The  equilibrium  proflle  is  expressed  os: 


where  A is  the  water  depth,  x is  the  horizontal  distance  from  the  ^ore.  and  A 
for  each  profile.  The  Bruun  rule  can  be  staled  in  a modifled  form  (Hands,  1981) : 
,_SWG, 


where  R is  shoreline  recession,  S is  sea  level  rise,  ir  is  the  width  of  the  active  portion  of 
the  profile  subject  to  adjustment,  G.  is  the  overfill  nuio  factor,  and  h,  is  the  vertical 
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distance  over  which  the  adjustment  occurs.  A schematic  representation  of  Bruun's  ruJe  is 

presented  in  ^£0^4-3. 

3.  Sediment  budget  method  (Everts.  1983):  The  sediment  budget  analysis  is  a method 
based  on  quantification  of  sources  and  sinks.  The  details  of  the  source  and  sink 
approach  are  documented  in  U.5.  Army  Corps  of  Engineers  Shore  Protection 
Manual  published  in  1984. 

4.  Dynamic  equilibrium  model  (Dean.  1983):  The  dynamic  equilibrium  model 
considers  changes  in  the  forcing  function  and  the  resuiting  transient  response 
characteristics  of  a beach  profile. 


Sea  level  afternse 
ImUal  sealevel 


Figure  4-3.  Schematic  representation  for  the  Braun  rule. 

Sea  level  rise  at  a particular  location  has  a number  of  components.  These 
components  include  global  (eustaiic)  sea  level  rise,  local  subsidence,  any  spatial 
variability  from  the  spatial  average,  local  meteorological  changes,  and  changes  in  the 
frequency  of  extreme  events  (Church  et  al..  IPCC  2001a).  Local  projections  of  sea  level 
rise  ore  thus  obtained  by  combining  the  local  factors  with  the  eustatic  sea  level  rise 
projections.  One  such  simplified  approach  suggests  extrapolating  the  historical 
observations  and  adding  global  eustaiic  projections  (Titus  and  Narayanan.  1993).  The 
method  suggested  by  Titus  and  Narayanan  assumes  that  the  global  sea  level  rose  12  cm 
over  the  last  century,  and  that  the  net  subsidence  at  a particular  location  was  1.2  mirt^ 
less  than  Lhe  observed  rate  of  relative  sea  level  rise  measured  by  ddal  gauges.  As  a result. 
Ihe  relative  sea  level  can  be  esiimaled  as: 


ioeol(r)  = gfoba/(r)+(r«W-  0.12)(/-1990) 


(4-3) 


In  equation  4-3. /ocflW  Is  the  rise  in  sea  level  by  year  t at  a poitieiilor  location 
(cm).  ghhalO)  is  the  global  sea  level  rise  projected  by  a particular  scenario,  and  trend  is 
the  current  relative  sea  level  rise  trend  at  a particular  location.  The  translation  ut 
shoreline  is  computed  in  terms  of  the  rise  in  reladve  sea  level  and  the  slope  of  the  coast. 
Schematic  representation  of  the  translation  is  presented  in  Figure  4-4. 


Figure  4-4.  Schematic  representation  of  translation  of  shoreline. 

The  translation  distance  T is  defined  as  (Hennccke  and  Cowell.  2000); 


where  dSL  represents  the  change  In  sea  level,  a represents  the  slope  of  the  coast,  and  T 
represents  the  translation  of  the  shore.  In  order  to  implement  Equation  4-4  in  grid 
analysis,  the  total  translation  at  any  transect  has  to  be  discretized  into  individual  cell  size 
uanslations  and  summed  up  for  a given  transect  to  obtain  the  total  translation.  The  total 
translation  along  a transect  is  thus  calculated  as: 

T = '^&T,  (4-5) 

where  jr,  represents  the  incremental  translation  at  cell  i of  a grid  with  cell  sizer:.  A 
schematic  representation  of  the  grid  aitalysis  of  translation  is  presented  in  Figure  4-5. 
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Figure  4-S.  Discrerization  ofincrrmenDl  translaiion  in  grid  anai^sis. 

The  siope  is  caicuialed  as  Ihe  maximum  rate  of  change  of  each  coll  from  Hs 
nei^bcra  using  the  Spatial  Artalysl  extension  of  ArcView,  The  algorithm  as  defined  by 

translation  for  cells  for  which  the  slope  value  is  zero  or  an  extremely  small  value  close  to 
zero.  In  order  to  avoid  the  numerical  error,  the  slope  of  a transect  is  calculated  as  the 
average  of  all  cells  along  the  trartsect.  This  method  ensures  numerical  stability  and 
provides  an  averaged  transect  slope  suitable  for  calculating  the  translation. 

During  n close  examination  of  the  methods  used  to  determine  translation,  it  became 
apparent  that  If  the  spatial  discmizailon  weresulficiaitly  small,  then  the  translation  fore 
given  cell  could  be  estimated  as  the  width  of  the  cell  along  the  transect.  The  final  value 
of  total  translaiion  along  a transect  could  be  computed  as  the  sum  of  contiguous  cell 
translations.  In  order  to  implement  this  method  one  has  to  code  the  algorithm  to  check 
for  each  cell  and  keep  track  of  Ihe  cells  ihot  meet  ihc  criteria  nnd  then  sum  their  values 


:t.  Therefore.  Ihc  cell  width  method  is  recommended  as  a more  stable 


method,  given 


sairsfred. 


IPCC  Climau  Change  Me 


rs 

The  most  recent  assessment  of  global  climate  change  and  its  impacts  conducted  by 
the  IPCC  was  released  in  200 1 . The  assessment  was  published  in  three  volumes  that 
summorisc  the  findings  of  the  three  working  groups  of  scientists.  The  contribution  of 
working  group  1 was  published  under  the  title  Climate  Change  2001;  The  Scientific 
Basis.  This  report  is  based  on  the  IPCC  Special  Report  on  Emission  Scenarios,  SRES 

greenhouse  gas  emissions  on  the  basis  of  a very  diverse  set  of  factors. 

The  scenarios  were  primarily  divided  into  four  families  (sioiylines).  I.e..  Al.  A2. 
Bl.  and  B2.  The  first  family  Al  is  further  divided  into  three  groups,  i.e.,  AIT,  AIFI,  and 
AIB.  All  SRES  scenarios  were  designed  as  quantitative  interpretations  (quantifications) 
of  the  SRES  qualitative  storylines.  Each  scenario  is  a particular  quantification  of  one  of 
the  four  storylines.  The  quantitative  inputs  for  each  scenario  involved,  for  instance. 

availability  of  various  forms  of  energy  agricultural  productivity,  and  local  pollution 
controls.  The  following  brief  description  is  reproduced  from  SRES  2000. 

The  A I storyline  and  scenario  family  describes  a future  world  of  very  rapid 
economic  growth,  low  population  growth,  and  the  rapid  introduction  of  new  and  more 
efficient  technologies.  Major  underlying  themes  are  conveigence  among  regions, 

reduction  in  regional  differences  in  per  capita  income,  "nte  Al  scenario  family  develops 
into  three  groups  that  describe  alternative  directions  of  technological  change  in  the 
energy  system 


V These  groups  are  AlFl,  AIB,  and  AIT. 


The  A2  storyline  and  scenario  family  describes  a very  heterogeneous  world.  The 
underlying  theme  is  self-reliance  and  preservation  of  local  identities.  Fertility  patterns 
across  regions  converge  very  slowly,  which  results  in  high  population  growth.  Economic 
development  is  primarily  regionally  oriented  and  per  capita  economic  growth  and 
technological  change  ore  more  fragmented  and  slower  than  in  other  stoiylincs. 

The  Bl  storyline  and  scenario  family  describes  a convergent  world  with  the  same 
low  population  growth  as  in  the  Al  storyline,  but  with  rapid  changes  in  economic 
structures  toward  a service  and  Information  economy,  with  reductions  in  material 
intensity,  and  the  introduction  of  dean  and  resource-efficienl  technologies.  The 
emphasis  is  on  global  solutions  to  economic,  social,  and  environmental  sustainability. 
Including  improved  e>Tuity.  but  without  additional  climate  inidalives. 

The  B2  sloiyline  and  scenario  family  describes  a world  in  which  the  emphasis  is  on 
local  solutions  to  economic,  social,  and  environmentiii  sustainability,  it  is  a world  with 
moderate  population  growth,  intermediate  levels  of  economic  development,  and  less 
rapid  and  more  diverse  technological  change  than  in  the  Bl  and  Al  stoiylincs.  While  the 
scenario  is  also  oriented  toward  environmemal  protection  and  social  equity,  il  focuses  on 
local  and  regional  levels. 

All  IPCC  scenarios  are  in  agreement  in  prediedng  an  increase  in  the  global  average 
icmpcrature  and  sea  level.  For  the  full  range  of  IFCC  scenarios,  the  global  sea  level  is 
projected  to  rise  by  0.09  to  0.S8  m between  1990  and  2100  with  a central  value  of  0.48  m. 
The  central  value  of  this  range  gives  an  average  rate  that  is  2.2  to  4.4  times  the  rate 
observed  during  the  20th  cenluiy.  Figure  4-S  shows  the  ranges  for  the  IPCC  scenarios. 


Global  sea  level  rise  estlmaies  and  a local  trend  values  are  used  to  obtain  reletlve 


sea  level  rise  values  for  a given  study  area.  The  global  estimates  serve  as  a basis  for  the 
lelallve  sea  level  rise.  However,  the  local  conditions  of  the  study  area  dictate  the  most 
probable  eslimaie  of  lalative  sea  level  rise.  Availability  of  longAerm  tidal  gage  records 


and  coastal  geonioiphology  data  improve  the  analysis  of  the  coastal  processes  governing 


the  rates  of  erosion  and  deposition. 


Figure  4-6.  IPCC  2001  Giobal  sea  level  rise  scenarios. 

Data  Availability  and  Qnality  Issues 

The  data  required  to  run  simulations  using  GENOEST  are  of  two  types.  The  first  is 
GIS  data  in  the  fonn  of  a DEM.  The  DEM  of  an  area  has  two  primary  parameters  that 
determine  the  precision  of  the  analysis.  These  are  the  vertical  resolution  and  the  cell  size. 
For  sea  level  rise  analysis,  the  vertical  resolution  of  the  DEM  is  preferred  to  be  floating 
pointdatawithabouta  meter  resolution.  The  cell  sin  of  the  DEM  detennincs  how 
precisely  the  horizontal  translation  and  area  submerged  are  to  be  calculated.  Most  DEMs 
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have  a cell  size  of  abou(  30  meters  or  more,  but  a cdl  size  of  about  10  meters  gives  a 
more  precise  estimate.  The  primary  problem  with  DEM  data  ofhighcr  vertical  and 
horizontal  resolution  is  the  size  of  the  data  file.  By  converting  a 30<meier  cell  size  data  to 
a 1 0-meter  cell  size  data,  the  file  dze  can  easily  triple.  The  la^e  data  size  requires  a 

[PCC  scenario  Is  readily  available  on  a decadal  basis.  If  the  analysis  requires  simulations 
for  smaller  time  steps,  then  the  data  series  can  be  obtained  by  consulting  the  IPCC 
headquaners.  The  second  approach  Involves  fitting  a curve  and  interpolating  or 
eittrapolaiine  as  required.  From  climate  modeling  point  of  view,  extrapolation  is  a not  a 
favorable  method  since  the  extrapolated  figures  are  based  purely  on  a mathematical 

estimates.  Attempts  to  extrapolate  the  data  for  another  100  years  have  resulted  in  a 
polynomial  set  of  equations  with  high  coefTicients  of  detemiinacion. 

The  polynomial  sets  of  equations  obtained  by  fitting  to  the  data  IPCC  forecast  data 
set  for  the  period  1990  to  2001  were  used  to  extrapolate  global  sea  level  rise  values  for 
the  period  2100  to  2200.  The  results  ofihe  extrapolation  and  the  corresponding 
coefficients  ofdeterminalion  are  presented  Figure  4-6  and  Table  4-1.  The  goodness  of  fit 
expressed  in  terms  of  coefflcieni  ofdeterminalion  is  well  above  (be  acceptable  limiL  In 
the  absence  of  a reliable  GCM  model  to  forecast  for  the  next  two  hundred  years  (he 
forecast  already  in  use  for  the  one  hundred  years  could  be  extrapolated  and  used  with 
reasonable  degree  of  acceptability. 
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Toble  4-1.  Goodness  of  til  for  polynomisl  curve  fining  of  the  )PCC  global  sea  level  rise 


CHAPTER  5 

DEVELOPMENT  OF  A VARIABLE  DENSITY  GROUNDWATER  ROW  MODEL 
TO  INVESTIGATE  THE  IMPACTS  OF  SEA  LEVEL  RISE  ON  COASTAL 
AQUIFERS 

BouDdorj'  CondilioDS  in  Croundwotcr  Flow  Modeling 
The  dcvelopmem  ofa  mimericol  groundwaler  flow  model  is  highly  dependenion 
Ihe  accurnlc  represcniation  of  Ihc  sysicm  boundary  condilions  (Andersen  and  Woessner, 
2002).  Careful  concepnializalion  of  bydrogeologic  boundaries  is  crilical  in  properly 
representing  the  system  using  mathematical  formuladons  (Wang  and  Anderson,  1982). 

system  boundary  condition  (Reilly,  2001 ).  Therefore,  the  modeler  has  to  methodically 
define  the  groundwater  flow  system  using  appropriate  governing  equations  and  identify 

The  equations  governing  groundwaler  flow  are  essentially  boundary  value 
problems  where  initial  and  boundary  conditions  need  to  be  specified  in  order  to  obtain 
solutions.  The  general  category  of  boundaries  in  flow  models  could  be  ulassilled  using 
the  common  mathematical  boundary  condilions  as  follows  (Reilly,  2001): 

1.  Dirichlei  boundary. 

2.  Neumann  boundary. 

3.  Cauchy  boundary. 

The  Dirichict  boundary,  also  known  as  Type  I boundary  or  specified  head 
boundary,  is  a calegoty  that  includes  specified  head  and  constant  head  boundaries.  The 
specified  head  boundary  is  a more  general  type  ofboundary  of  which  the  constant  head 
boundary  is  a special  case.  The  specified  head  boundary  requires  Ihe  hydraulic  head  h to 

78 


79 

(he  head  U held  com(am  al  a givea  location  in  space  and  time  (Franke  et  al,  1 987).  In 
reaUly.  a conaiani  head  boundaiy  (also  rereired  to  as  constanl-potemial  or  eijuipoletiuai 
boundary  I occurs  where  a pan  or  the  boundary  surface  of  an  aquifer  system  coincides 
with  a surface  of  spatially  and  temporally  urriform  head  values. 

The  Neumann  Boundary  also  known  as  a Type  2 boundary  or  specified  (low 
boundary,  is  a category  that  Includes  specified  flux  and  no-flow  boundaries  (also  known 

where  flux  q (discharge  per  unit  cross-sectional  area)  varies  with  space  and  time.  The 
flux  could  also  be  sfiecified  os  constant  in  time  but  variable  in  space.  A no-flow 
boundary  on  (he  other  hand  is  a special  case  of  the  specified  flux  boundary  where  (he  flux 

The  Cauchy  Boundary.  aJso  known  as  Type  3 or  head  dependent  flow  boundary, 
includes  the  head  dependent  flux  boundary.  This  type  of  boundary  is  encountered  when 
the  flow  in  an  aquifer  changes  in  response  to  the  variability  in  head  of  the  adjacent 

A free  surface  boundary  is  characterized  by  a change  in  the  head  value  with  respect 
to  elevation  through  time.  i.e_  fr  = z or,  mote  generally,  h = /z).  This  type  of  boundary  is 
used  to  represent  a water  table  and  the  transition  between  freshwater  and  underlying 
seawater  inacDOSlalaquiferfFrankeetal.,  1987).  Ifthediffusive-dispeisive  process  of 
solute  transport  is  neglected  and  the  saline  groundwater  located  in  Ihe  intruded  seaward 
zone  Is  assumed  10  be  static,  then  ihe  frcshwaier-sollwaler  iransiiion  zone  can  be  treated 
as  a sharp  interface.  The  assumption  of  a shoip  inicrface  allows  (be  interface  to  be 


rfflce  (no-flc 


ndary) of the 


i groundwate 


flow  system  (Franke  el  al„  1987). 

A seepage  surface,  or  seepage  face  boundary,  is  encountered  when  a saturated  flow 
mathematical  representation,  the  seepage  face  issimiiarto  the  free  surface  where  the 
A genertdized  classification  of  the  major  boundaiy  types  (after  Franke  et  ol..  1987 


and  Reilly,  2001)  is  presented  in  table  M . 

Table  5-1  Ciassifleation  of  boundary  condii 
Boundary  Type  Boundary  Name 


Type  1 Specified  head 

Constant  head 
Type  II  Specified  flux 

Stream  line 

Type  III  ITead  dependent 

flux 

Free  surface 
Seepage  face 


Name 

Dirichlel 

Neumann 


Long-Term  Response  of  Coastal  Groundwater  Flow  to  Tidal  Fluctuations 
Coastal  aquifers  are  usually  chamcierured  by  the  dynamic  interaction  between 
seawater  and  groundwater,  which  is  influenced  to  a certain  degree  by  the  fluctuation  of 

coosiaJ  aquifers  have  shown  a lagged  fluctuation  in  piczometric  levels.  In  almost  all  of 
the  studies,  this  response  is  observed  to  be  signiricanl  on  a smaller  temporal  scale, 
primarily  on  an  hourly  basis  (Erskine.  1991:  White  and  Roberts,  1994;  Millhamand 
en  and  Jiao.  1999;  Jiao  and 


Howes,  I99S;  Che 


I Tang,  1999).  Tidal  infill 


81 

graundwam  level  is  also  dependenl  on  the  proximity  of  the  observation  point  to  the  shore 
(Jiao  and  Tang,  )999).  As  pan  of  the  research  conducted  for  this  dissertation,  studies 
were  conducted  using  tidal  and  aquifer  data  for  south  Florida.  Theaimofihc  research 
was  to  determine  whether  tidal  data  should  be  part  of  the  boundary  condition  for  the 
groundwater  flow  model. 

The  response  of  coastal  aquifers  to  tidal  influence  also  depends  on  the  type  of 
aquifer.  There  is  a marked  difference  in  response  to  tidal  fluctuation  between  confined 
and  unconflned  aquifers  [White  and  Roberts,  1994).  The  spatial  extent  of  the  response 
aone  measured  from  the  open  water  may  vary  from  a hundred  meters  or  more  for 
confined  aquifers  to  about  30  or  40  meters  for  unconfmed  aquifers  (Jiao  and  Tang,  1999). 
Although  unconfmed  aquifers  exhibit  higher  attenuation  and  negligible  response,  the 
leakage  from  overlying  unconTined  aquifers  may  affect  die  response  of  underlying 
confined  aquifers.  For  the  purpose  of  illustrating  the  amplitude  of  tidal  fluctuation  and 
groundwater  head  in  coastal  oquifeis  a schematic  cross-section  ofa  coastal  aquifer  with 
tidal  cfTcci  on  the  polentlomelric  surface  is  presented  in  Figure  S-l. 


Figure  5-1.  Coastal  aquifer  with  tidal  effect  on  the  poiemiomelric  surface. 

The  tide  varies  between  two  extremes  of  high  and  low.  The  time  elapsed  for  the 
tide  to  vary  between  these  two  extremes  is  referred  to  as  tidal  period.  For  a tide  with 


amplitude  of  tidal  change  and  a tidal  period  the  amplitude  of  tidal  nuctualion  Is 
given  aa  (Jacob,  I9S0): 


(5-1) 


where  Ht  is  the  ampUtude  of  tidal  fluctuation,  a parameter  that  represents  the  variation  in 
the  piezometric  head  observed  at  the  perpendicular  distance  x off  the  coast  due  to  the  tide 
at  any  distance  x measured  inland  from  the  coast.  The  amplitude  of  tidal  fluctuation  was 
plotted  for  varying  aquifer  mnstnissivity  values  and  presented  in  Figure  5-2  using 
storage  of  0.001  and  tidal  period  of  0,25  day  (d  hours). 


i 

I 


Figure  5-2,  Amplitude  of  tidal  Ouctuaiitm  in  a confined  aquifer. 

The  response  of  coastal  aquifers  to  tidal  fluctuatimtsis  not  Instantaneous,  but 
rather,  it  is  lagged  by  some  time.  Referring  to  Figure  5-1,  for  an  aquifer  ofstorativity  S 
and  transmissivity  T.  the  time  lag  t.  is  given  by  (Fetter,  2001): 


From  Equations  5-1  and  5-2  it  is  clear  that  the  response  of  aquifers  to  tidal 
(luclualion  exhibits  exponential  decrease  as  the  location  gets  farther  away  from  the  coast. 
Similarly,  the  time  lag  decreases  linearly  with  distance  from  the  coast.  Also,  as 
demonstrated  in  Figure  5-2,llte  amplitude  of  the  tidal  fluctuation  increases  with  increase 
in  transmissivity  of  the  aquifer. 

The  form  of  the  relationship  that  exists  between  piezometer  readings  and  tidal 
fluctuations  is  approximated  to  be  sinusoidal.  Once  the  lime  lag  between  piezomelric 
level  and  tidal  fluctuation  is  computed,  then  it  is  possible  to  superimpose  the  plot  of  one 
on  another  by  applying  a time  lag  and  a tidal  efficiency  factor  (Erskine,  199i|.  Fora 
homogeneous  conRned  aquifer  sinusoidal  oscillotions  bt  pressures  propagate  along  the 
aquifer  according  to  the  following  equation  (Ferris,  1951): 

<S-3) 

where  A = the  groundwater  head  relative  to  mean  sea  level,  x = distance  from  sea.  < ■ 
lime.  r,  = period  of  tidal  oscillation,  A,  = amplitude  of  tidal  oscillation,  transmissivity 
ofaquifer.  and  5 = aquifer  slomtiviiy.  The  time  lag  is  defined  as  explained  in  Equation 
(5-2)  and  the  tidal  efficiency  factor  t;  is  defined  as  (Erskine,  1991); 
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Meaascalr^el  ■ 


Mt«n  frouiulw^  level' 


UncouBped  aquifer 


Leaky  cooEiwd  aquifer 


Figure  5-3.  Leaky  confined  aquifer  with  tidal  Influence. 

Jiao  and  Tang  (1999)  derived  on  equmion  governing  the  flucluaiion  in  hydraulic 
head  due  to  tidal  fluctuation  for  a leaky,  confined,  hotnogeneoua,  and  isotropic  aquifer 


from  a specified  datum  as  shown  in  Figure  5-3  (Jiao  and  Tang.  1999).  With  the 
assumption  of  a negligible  storage  of  the  semi  permeable  layer,  the  lealcage  is 
proportional  to  the  difference  in  head  between  the  two  aquifers.  For  the  specified  set  of 
conditions  a one-dimensioiiol  gioundwaier  flow  equation  is  defined  as  (Jiao  and  Tang, 
1999): 


tlie  head  boundary  condition  is  defined  at  the  tidal  boundary  of  the  coast  as: 


governed  by  Dupuil's  assumption  with  a i 


‘ head  of  h.  measured 


(5-5) 


(5-6) 


Further  inland,  i 


Ji(0.r)  = 6,  +.4eos(iuM-c) 

CO , the  head  boundary  condition  is: 
6(co,r)  = A, 


(5-7) 


where  ^ Is  Ihe  amplitude  of  [he  lidaj  change,  c is  the  phase  shift,  and  w Is  the  tidal  speed 
defined  as  2s/lo.  With  the  boundaty  conditions  specified,  the  solution  to  equation  5-5  is 
(Jiao  and  Tang.  1999): 

h(i,i)eA,a-.4e)ip(-pa)coa^a-f-^+cj  (5-8) 

where  p is  a term  defined  as: 


for  an  aquifer  with  zero  ieakance  then  the  solution  reduces  to; 

(5-10) 

in  situations  where  the  proximity  of  the  wells  and  the  length  of  time  period  of 
simulation  wotranls  the  use  of  tidal  data,  corrections  need  to  be  applied  to  compensate  for 
tidal  crfecls.  Thecorrectionsareappliedbyfilteriogtbeobservedpiezometricdata.  In 
order  to  filter  the  data,  tidal  efficiency  and  time  lag  have  to  be  determined  using 
Equations  5-2  and  5-4. 

For  a number  of  observation  piezometers  located  at  distances  ranging  from  50  m to 
400  m from  the  seashore.  Etskine  (1991)  used  the  following  lillering  equation: 

A,.(0  = M0-i7[r(f-(,)-f]  (5-U) 

where  A,  (t)  • filtered  piezomeiric  level  at  lime  / (m);  A(i)  = piezometric  level  at  lime 


i(m);  f • mean  tide  level  (ra);  r(/)“tide  level  atlimet  (m). 


As  port  of  this  study  select  tidal  stations  and  monitoring  wells  were  analyzed  in  the 
Slate  of  Florida.  The  locations  and  names  of  the  stations  are  presented  in  Figure  5-4  and 
Figure  5-5.  The  tidal  stations  data  was  obtained  from  the  Center  for  Operational 
Oceanographic  Products  and  Services  (COOPS),  which  is  a data  dissemination  center  of 
the  National  Oceanic  and  Atmospheric  Administration  (NOAA).  The  groundwater  level 
data  was  obtained  Irom  the  US  Geological  Survey  and  South  Florida  Water  Managemeni 
DistricL 

From  the  available  tidal  stations  for  the  state  of  Florida,  those  in  the  area  of  interest, 
which  primarily  consists  of  three  counties  in  south  Florida,  were  selected  and 
corresponding  groiindwaler  monitoring  sites  were  selected  based  on  proximity.  Although 
it  was  difficult  to  find  corresponding  sites  where  both  a tidal  station  and  groundwater 
monitoring  site  were  located  close  to  each  other,  it  was  possible  lo  find  some  sites  where 
the  tidal  stations  and  groundwater  monlloring  slleswere  within  a few  thousand  meters  of 
each  other.  Shorter  time  period  lime  series  and  daily  average  tidal  data  for  the  year  2003 
were  plotted  for  the  selected  (see  Figure  5-6).  Corresponding  daily  groundwater  level 
data  also  were  plotted  for  the  corresponding  groundwater  sites  (see  Figure  5-7), 
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Figure  S-6.  'Hdal  data  of  select  stations  in  south  Florida.  A)  Boynton  Beach.  Culler 

Biscayne  Bay  and  Hobe  Sound.  B)  Key  West.  Miami  Beach,  North  Miami 
Beach.  C)  North  Palm  Beach.  Palm  Beach,  Vaca  Key. 
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Figure  5-«.Cc 


Figure  5-6.  Continued 


luiconnned  aquifer  In  the  area.  The  water  level  dau  of  monitoring  well  G-3327_G. 


flucluallon  was  computed  using  Equation  S-3  (Ferris,  l9SI)and  thecalibraiion  daiaof 
aquifer  parameters  (or  the  Biscayne  Bay  regional  model  (Langevin.JOOl).  For  the  sake 


ofcompanson.  the  groundwater  head  fluctuation  was  also  computed  again  using  the 
EguationS-IO  (Jiao  and  Tang.  1999).  The  data  for  both  cases  consisted  oTa  surficial 
aquifer  of  average  thickness  4S  m.  a hydraulic  conductivity  9000  m/day.  storage  of  0.001. 
and  an  average  tidal  fluctuation  of  with  varying  periods  of  tidal  oaclltallon  that  ranged 
from  0.5  hours  to  24  hours. 


The  results  iTom  Equation  S-3  showed  almost  no  fluctuation  even  at  the  coast,  as 
the  fluctuation  was  limited  to  less  than  SOm.  The  results  5om  Equation  5*10  showed 
small  fluctuations  up  to  a distance  of  200  m.  but  after  that,  the  oscillation  becomes 
insignificant.  At  a distance  of  12.000  m.  there  is  no  Influence  on  the  groundwater  head. 

The  calculations  made  to  investigate  the  influence  of  tidal  fluctuations 
demonstrated  that  theefibcts  tidal  fluctuations  have  on  groundwater  levels  over  a longer 
period  of  simulation  that  ranges  from  decades  to  centuries  could  be  neglected.  Instead, 
the  focus  should  be  on  the  sustained  rise  in  sea  level.  The  scale  of  a model  also  plays  a 
signiflcanl  role.  Groundwater  head  in  areas  beyond  a couple  of  hundred  meters  in  the 
case  of  confined  aquifers,  and  perhaps  less  than  fifty  meters  in  the  case  of  unconfined 
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squifotsdoes  noi  seem  to  be  influenced  by  tidaJ  flueiuntion.  Usually  regional  models 
have  spatial  exienls  of  a thousand  meters  or  more,  thus  the  influence  is  not  significant  to 
be  considered  due  to  ancnuoiion. 

The  spatial  and  temporal  scale  of  groundtsaier  flow  models  varies  depending  on  the 

requirements.  If  the  spatial  scale  of  the  models  is  very  small,  the  transect  perpendicular 

to  the  seaside  boundary  is  relatively  short,  and  the  timescale  of  modeling  is  also 
relatively  small,  then  consideration  of  tidal  influence  might  be  required.  In  situations 
where  tidal  influence  has  to  be  considered,  coireclioiis  need  to  be  applied  to  the  observed 
piezometric  readings.  Since  most  regional  aquifer  models  hove  cell  sizes  of  a few 
hundred  meters,  the  adjoining  cells  at  the  seaside  should  be  defined  as  a Dirichlet 
boQiidaiy.  and  the  head  and  concentration  should  be  specified  accordingly.  For  the  head 
boundary,  the  value  should  he  set  to  vaiy  iu  response  to  sea  level  rise. 

ModificaHon  of  Boundary  Conditions  to  Account  for  Impacts  of  Sea  Level  Rise 

The  hydrogeologic  boundaries  of  coastal  aquifers  are  complex  due  to  the  presence 
of  a seaside  boundary  in  addition  to  the  inland  boundaries.  During  the  investigation  of 
seawater  intrusion  and  sea  level  rise,  the  primary  boundary  condition  of  concern  is  the 
seaside  boundary.  The  more  customary  way  ofspecilying  the  seaside  boundary 
condition  involves  defining  a Dirichlet  boundary  condition  where  the  head  is  specified  or 
held  constant.  For  a variable  density  model,  concemnuion  values  at  the  seaside  are  also 
held  consuuu.  In  order  to  define  the  change  in  the  head  boundary  during  transient 
simulations,  numerical  codes  utilize  a number  of  techniques.  For  example,  in  the  USGS 
code  MODFLOW.  a linear  relationship  can  be  assumed  between  specified  starting  and 
final  heads  for  each  stress  period.  In  between  stress  periods,  the  specified  head  boundaiy 
values  are  intctpolated  linearly  between  the  specified  starting  and  ending  head  values  for 
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Ihe  different  time  steps  within  a given  stress  period.  During  a given  time  step,  the  head 
values  are  held  constant  and  the  interpolation  algorithm  is  executed  once  every  lime  siep. 
SEAWATSLR:  ModincalioB  of  SEAWAT  to  Model  Sea  Level  Rise  Impacts 
The  now  part  of  the  three-dimensional  variable  density  groundwater  now  modeling 
code  SEAWAT  is  based  on  MODFLOW.  As  described  in  Chapter  2.  SEAWAT  is 
formulated  to  yield  values  of  groundwater  head  In  terras  of  equivalent  freshwater  head. 
One  of  the  How  packages  incorporated  in  SEAWAT  is  the  rime  Variant  Specified  Head 
package  (CHD).  The  CHD  package  is  used  to  represent  hydrogeologic/hydrologic 
features  with  static  or  nearly  static  head  values.  Customarily,  Ihe  seaside  boundary  is 
incorporated  into  a groundwater  flow  model  as  a conslam  head  boundary.  The  water 
hydraulic  head  on  the  seaside  is  expressed  as  equivalent  freshwater  head  with  hydrostatic 

The  modified  version  of  SEAWAT  named  SEAWATSLR  accounts  for  lime  variant 
changes  in  the  seaside  boundary  by  incorporating  incremental  values  of  rise  in  sea  level 

for  each  flow  time  step.  The  changes  in  lime  variant  constant  head  are  thus  expressed  as 

boundary  changes  due  to  suainined  sea  level  rise  through  modification  of  the  CHD 
algorithms. 

SEAWATSLR  incorporates  a modified  linear  inleipolalion  routine  in  the  CHD 

sea  level  specified  by  the  user.  The  incremental  sea  level  value  is  obtained  from  Ihe 
relative  sea  level  rise  modeling  results  obtained  by  using  CCGEN,  as  described  in 
Chapter  4.  The  sea  level  rise  estimated  for  a given  forecast  period  is  distributed  as  an 


96 

incremcnial  increase  and  included  in  the  conslanl  head  for  Ihe  entire  groundwater  flow 
simulation  period  at  every  time  step. 

Tlie  etiuivalentfreshwaterhead  in  thelime  variant  specified  head  boundary  is 
written  for  every  cell  (y,*)  at  any  time  I as: 


Pf 


(S-I2I 


Transient  changes  in  ihe  boundary  head  conditions  result  in  a variation  in  the  head 
term  between  successive  time  steps.  The  head  at  ihe  nei«  time  step  is  related  to  Ihe  head 
at  the  current  time  step  as  follows: 


*"Vo.w)  =*'/(,. ,u,  + iv  dr  (5-13) 

where  (V  is  the  rate  of  sea  level  rise  and  dr  is  the  elapsed  time  at  the  given  flow  lime 
step.  Substituting  for  the  head  tern  in  lime  step  1+ 1,  the  etjuivalent  freshwater  head  for 
Ihe  lime  variant  conslanl  head  boundary  becomes 


The  head  between  two  time  steps  is  interpolated  linearly  in  the  current  formulation 
of  SEA  WAT.  However,  the  seaside  boundary  which  in  coastal  aquifers  is  related  to 
sustained  sea  level  has  a non  linear  variation  through  time.  Referring  to  the  best  fit 
equations  developed  to  extrapolate  the  global  sea  level  rise  scenarios  in  chapter  4.  the 
relationship  has  been  determined  to  be  polynomial.  According  to  the  current  formulation 
of  SEA  WAT.  the  head  is  linearly  interpolated  as 


(5-15) 


where  4/,  is  a time  fraction  expressed  as  the  proportion  of  stress  period  to  the  end  of  the 

current  time  step.  This  lime  fraction  iscalcuiatcd  as  the  ratio  of  total  ciapsed  time  to  the 
iength  of  stress  period 


PERTIM 

PERLEN 


(5-16) 


where  PERTIM  is  the  totai  time  elapsed  during  a stress  period  and  PERLEN  is  the  iength 
of  the  stress  period.  With  the  current  definition,  equation  5-14  is  updated  to  utiiize  the 
accumuiated  time  in  a given  time  step  4/  = PERTIM  mA'A  becomes 


The  time  component  of  the  term  that  accounts  for  the  rale  of  sea  ievelriseis 
customariiy  expressed  as  per  century  or  per  year.  Therefore,  the  user  has  to  accordingly 
change  it  to  days  which  is  the  more  common  time  scale  used  in  groundwater  mode  iing. 
AJso  for  a century  iong  simulation  using  a time  scale  less  than  a day  will  make  the 
computational  efTorts  unnecessarily  high.  The  length  of  a time  step  in  a stress  period  Is 
defined  by  the  number  of  time  steps  and  the  time  step  multiplier,  and  the  length  of  time 
step  is  defined  as: 

PERLEN(.\-TSMULT) 

H-TSMULT)'^ 

where  DELT  is  the  time  step  length,  TSMULT  is  the  lime  step  multiplier,  and  NSTP  is 
the  number  of  time  steps.  Using  the  concept  of  piecewise  linearization,  the  linear 
interpolation  scheme  implemented  in  SEA  WAT  can  be  used  without  modification.  The 
piecewise  linearization  within  the  specified  time  slep  of  simulation  is  an  acceptable 
approximation  within  the  specified  time  slep.  which  is  usually  smaller  than  the  century 


flowchan  for  SEA  WATSLR  (modified  Irom  SEA  WAT)  is  prcsenlcd  In  Figure  5-9. 
Procedures  in  red  boxes  show  components  modified  or  added  o.s  part  of  this  research. 


>rSEAWATSLR. 


CHAPTER6 

DEVELOPMENT  OF  EVOLUTIONARY  ALGORITHM  BASED  GROUNDWATER 
FLOW  OPTIMIZATION  MODEL 

FormuloliOD  of  Objective  Function  and  Constraints 
The  objective  functions  to  be  solved  in  oplimizalion  ptobleius  of  coastal  aquifers 
usually  involve  maximizing  the  benefits,  minimizing  the  cost  of  pumping,  minimizing  the 
advancement  of  saltwater  farther  inland,  minimizing  drawdown  due  to  pumping  below  a 
specified  value,  and  maximizing  the  pumpage  (Shamir  el  al„  1984;  Willis  and  Finney, 
1988;  Finney  at  al„  1992;  Hallaji  and  YazicigiL  1 996;  Emch  and  Yeh,  1998;  Das  and 
DaRa,  1999a,  1999b;  Cheng  et  al„  2000;  Benachmi  el  al.,  2000).  Most  optimization 
approaches  used  in  coastal  aquifers  management  focus  on  maximizing  pumpage  and 
minimizing  cost  subject  to  the  eonstrainu  of  drawdowns  and  saltwater  intrusion  criteria. 

In  this  research,  two  approaches  were  considered  for  the  purpose  of  demonstration, 
hut  only  one  was  implemented  in  the  development  of  the  optimization  model.  The  first 
approach  is  a sharp  interface  abroach  (Emch  and  Yeh,  1998).  while  the  second  iqipioach 
is  a variable  density  approach  (Gordon  et  al.,  2001).  The  sharp  interface  approach  was 
modified  to  work  with  the  formulation  ofSEAWATjusl  for  the  sake  of  demonstrating 
the  feasibility  of  the  approach  for  future  research.  The  variable  density  approach  was 
applied  using  genetic  algorithms  and  the  results  arc  presented  in  Chapter  7.  Both 
approaches  are  based  on  a multi-management  period,  mulb -objective  optimization 
problem  formulation,  but  each  one  ubiizes  a different  approach  of  constraim  and 
objective  function  formulation. 
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Muili'Objeclive  optimization  problems  require  finding  a vector  of  decision 
variables  that  satisfy  constraints  and  optimize  a vector  function  whose  elements  represent 
the  objective  function  {Oscycka,  1935).  For  a vector  of  rt  decision  variables  denoted 
mathematically  denoted  asx,,i  ■ l,n,  the  vector  noution  of  decision  vanables  is  defined 
as  (Cocllo  Coello  et  al.,  2002): 


The  constraints  that  need  to  be  satisfied  are  defined  are  expressed  as  inequalities: 


where  the  number  of  equality  constraints  p is  less  than  the  number  of  decision  variables  n 
in  order  to  avoid  the  formulation  of  an  over'consirained  problem.  The  number  of  degrees 
of  freedom  is  then  defined  as  fr-p,  thus  ensuring  the  number  of  unknowns  is  less  than  the 
number  of  equations. 

The  objective  functions  are  defined  afi(i)j2(x)....f^x)  where  k is  the  number 
of  objective  functions  in  the  multUobjcctivc  problem  being  solved.  The  objective 
fiuiclions  arc  represented  io  vector  form  os: 


(6-1) 


(6-2) 


or  as  equalities  as: 


(6-3) 


The  general  multi' 


I is  defined  in  Euclidean  n-space 


(6-5) 


which  satisfies  the  m inequality 

<=l <".»,[x]  = 0 t = l p 


and  uplimizes  the  vector  function: 


?(1]. 


<=)1 

r,p] 


(6-6) 


Sharp  Interface  Approach 

The  first  approach  is  based  on  the  position  of  the  saltwater  freshwater  interface  and 
fonnulates  the  objective  function  on  the  basis  of  minimizing  the  cost  and  volume  of 
saltwater  in  the  aquifer.  Theotigmal  formulalion  (Entch  and  Yeh.  1998)  used  a sharp 
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inierface  approach,  and  ihe  flow  model  wai  simulated  using  the  USGS  code  SHARP 
(Essnid,  1 990a).  The  oplimizalion  was  performed  using  the  commercial  software 

MINOS  (Muitagh  and  Saunders.  19g7).  In  order  to  solve  the  nonlinear  objectives  and 
constraints,  a projected  augmented  Eulerian  Lagrangian  algorithm  wa.s  used  and  the 
nonlinear  problem  was  reduced  to  a series  of  linearly  constrained  sub-problems.  The 
reduced  sub-problems  were  then  solved  using  the  reduced  grtwlienl  algorithm. 

The  formulation  ofEmch  and  Yeh  was  modified  for  this  discussion  and 
reformulated  on  Ihe  basis  of  cquivaJem  freshwater  heads,  the  position  of  the  saltwater- 
freshwater  interface,  and  freshwater  withdrawal  through  pumping.  The  simulation  and 
opumizalion  models  could  be  coupled  using  the  response  matri.s  approach.  This 
approach  requires  generating  the  unit  responses  of  aquifer  to  pumping  stresses  and 
utilizing  these  responses  in  the  objective  function  lo  be  optimized. 

Minimizing  the  cosi  of  pumping  and  maximizing  the  amount  of  water  pumped 
require  multi-objective  crileria,  thus  satisfying  the  condition  for  a muili-objective 
optimization  model.  The  objective  function  could  be  formulated  as  a multi-objective, 
mulli-management  period  optimization  problem. 


Figure  6-i.  Represeniaiion  of  a coascal  aquifer  for  objective  function  formulation. 
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Consider  Ihe  single  layer  aquifn  shown  in  Figure  6<1.  Let  G be  the  set  of  all 
layered  aquifers  front  which  any  layer  i ! G is  considered.  From  conrinirity  of  fluid 
pressure  the  elevation  oflhe  saltwater-freshwater  interface  above  the  aquifer  datum  is: 

<6-7) 

where  and  A'  are  vertically  averaged  freshwater  and  saltwater  heads;  pj  and  p^  arc 
freshwater  and  saltwater  densities;  p'^  and  p\  are  freshwater  and  saltwater  density 

ratios;  $- — — — .and  Z,  is  the  elevation  ofinterface  above  datum. 

P.-Pl 

Assume  there  are  n decision  variables  in  the  problem  which  are  represented  in 
vector  foim  as: 

F = (/>, P.,) 

where  at  = at,  +ct,  is  the  total  number  of  supply  sources  and  n s:  a>  k r . 

The  state  variables  in  the  formulation  suggested  are  freshwater  heads,  saltwater 
heads,  and  the  position  of  the  interface.  In  order  to  apply  this  approach  using  SEA  WAT, 
the  equivalent  freshwater  head  and  the  position  of  the  interface  could  be  used.  This 
requires  isolating  a specified  isochlor  as  the  defining  interface  boundary  and  based  on 
that  calculating  the  location  of  the  Z,  boundary. 

The  cost  objective  requires  minimising  the  operation  cost  of  pumping  not  including 
the  cost  of  drilling  the  wells  and  constructing  the  supply  system.  The  cost  minimization 
could  be  formulated  as: 


(6-9) 
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where  Z,  “elevation  of  the  interface,  P,_,  “ water  supply  rale  from  souree  i for  time 
period  j.  C,  ■ unit  cos!  of  groundwale  extraction  per  heighl  of  required  lift  or  of  surface 
water  supply  for  source  t.  ft,  = set  of  groundwater  sources,  and  Q,  is  set  of  surface 
water  sources 

The  saltwater  intrusion  objective  function  was  formulated  in  terms  of  the  volume  of 
saltwater  in  the  aquifer 


The  output  ftom  SEAWAT  lain  terms  of  the  equivalent  freshwater  head  therefore 
Equation  (6-9)  was  refortimlalcd  accordingly.  The  position  of  the  interface  is  similar  to 
the  saltwater  head  representation  in  SEAWAT.  thus  Z\  could  be  replaced  by  A,,  which 
was  converted  from  the  equivalent  freshwater  head  formulation- 
The  reformulated  objective  ftinction; 


Mn  Z,  = Jljjz;(.c,y)<ira'y| 


(6-10) 


Min  Z,  =|][j|A,^(x,y)iirfl'yJ 


(6-11) 


substituting  for  A,', , 


(6-12) 


and  the  constraints  include: 


Demand  constraint: 


(6-13) 


Well  capacity  constraint: 


Ose,^Sg,*“  Vied, 


(6-14) 
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Drawdown  constiaints: 

(6.15) 

Surface-water  supply  consnainls: 

Ose,^se“  Vien,  (6-16) 

Variable  Density  Approacb 

The  second  approach  is  based  on  the  forniulation  of  a regional  scale  opiimiaalion 
model  (Gordon  el  a].,  2001).  The  objeclive  oflhis  optimizalion  mode!  is  maximirini  Ihe 
amounl  of  water  pumped  and  minimizing  the  amouni  of  sail  mass  ezuacicd  with  the 
water.  Thia  requires  performing  a mass  balance  for  Ihe  flow  simulalion  and  calculating 
(he  amounl  of  sail  mass  exuncted  for  each  managemenl  period.  The  response  matrix  was 
gcneraled  from  pumping  sliesses,  mass  of  sail  extracted  through  pumping,  and  the 
constraints  on  pumping  for  each  managemenl  period. 

Consider  a regional  aquifer  with  a number  ofpumping  wells,  ll  is  assumed  thol  the 
aquifer  is  already  subject  to  saltwater  intrusion  and  the  objective  is  pumping  the 
maximum  amount  of  water  without  violating  a set  of  predefined  waler  quality  and 
quantity  consuainls.  Let  Ni,  be  Ihe  total  number  of  wells  and  N,,  be  Ihe  number  of  lime 
periods  in  the  entire  managemenl  period  during  which  the  pumping  stresses  are  constant 
In  the  current  formulation  these  time  periods  are  set  to  be  Ihe  same  as  the  stress  periods 
of  the  flow  model.  The  decision  variables  of  the  management  model  are  the  pumping 
rales  with  consuaims  on  well  pumping  rales  and  total  amount  of  water  extracted. 

The  mathematical  representation  of  the  governing  ntanagemeni  model  is 
formulated  as  follows  (Gordon  el  al.,  2(X)I); 
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Max 


(6-17) 


(6-18) 


e_,se«so 


(6-19) 


(6-20) 


where  is  Ihe  pumping  rale  of  well  ip  in  time  period  ;/j  (in  this  case  stress  period)  and 

is  constant  and  considered  negative  to  be  consistent  with  the  formulation  of  groundwater 
models;  (r)is  the  concentration  of  salt  in  welli^  during  time  period  r;  is  the 
maximum  pumping  capacity  of  well  yt;  and  go„  and  go^  are  the  lower  and  upper 
bounds  of  the  total  amoiml  of  water  tbai  can  be  extracted  from  the  aquifer. 

The  presence  of  the  concentration  term  in  the  objective  function  makes  the  problem 
defined  above  non-linear,  non-convex.  and  non-smooth.  The  solution  adopted  by  Gordon 
el  al.,  (2001)  was  the  Bundle-Trust  method,  which  is  a modification  of  the  classical 
gradient  methods.  The  constraints  defined  in  (6-19)  are  incorporated  directly  into  the 
opiimixslion  model  while  the  consiiainls  defined  In  (6-20)  are  incorporated  as  e penally 
term  giving  a modified  management  model  with  a weighted  sum  of  the  two  objective 
functions.  The  modified  objective  function  thus  becomes: 


subject  to  S s 0 . The  weights  Pi  and  P;  are  adjusted  to  make  sure  that  the 
upper  and  lower  ptunping  bounds  are  not  violated. 

CENOPT:  A Genetic  Algorithm  Based  Solution  of  Groundwater  Flow  OpUmiaalion 
Method 

The  formulation  of  practical  water  resources  opiimiaatton  problems  is  often 
characterized  by  conveniot  and  non-linearity.  As  a result,  obtaining  second  order  and 
even  in  some  eases  first  order  derivatives  of  the  objective  functions,  which  are  required 
by  many  traditional  optimization  methods,  is  di  Bicult  if  not  impossible  (Tang  and  Mays, 
I99S).  Moreover,  it  is  difficult  to  ensure  that  the  solution  obtained  is  the  globally  optimal 
solution  for  the  given  problem.  As  eaplained  in  Chapter  3.  Genetic  Algorithms  differ 
from  the  traditional  optimization  methods  in  many  ways  and  are  capable  of  handling 
conveit  and  non-linear  problems  by  ensuring  global  or  near  global  optimal  solution. 

Once  the  objective  function,  the  decision  variables,  and  the  constraints  are 
formulated,  then  the  decision  variables  are  encoded  into  strings  referred  to  as 
chromosomes.  Commonly  these  strings  ore  coded  as  binary  strings,  but  it  is  also  possible 
to  encode  them  as  integer  or  real.  Once  the  encoding  is  accomplished,  then  a random 
number  generation  algorithm  is  used  to  initialize  the  first  population  of  possible 
solutions.  Equal  probability  of  survival  is  assigned  initially  for  the  entire  ensemble.  The 
objective  function  is  used  to  define  the  fibtess  function  of  the  problem.  The  fitness 
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function  assigns  a probability  for  each  chromosoms  to  generate  new  chromosomes  and 
then  tests  iu  performance  in  between  rwo  successive  generations  by  introducing 
selection,  crossover,  and  mutation. 

The  basic  features  of  the  core  GA  operations  include  creep  mutations,  elitism, 

niching,  tournament  selection,  single  point  or  imifomi  crossover.  The  detailed  operations 

of  the  GA  operation  are  presented  in  the  process  flowchart  figure  6-2,  the  processes 
generally  include  the  following: 

I . Randomly  generate  on  initial  seed  population  of  binary  strings  representing 

possible  designs.  This  initial  population  of  strings  represents  realiaations  of  a set  of 
decision  variables.  Each  string  is  coded  as  binary  so  that  rhe  alleles  (particular 
locations  on  the  strings)  will  have  values  of  Os  andls. 

2-  Calculate  the  fitness  of  each  string  of  the  current  generation  using  the  objective 
function  and  the  constraints,  The  strings  are  then  ranked  using  their  fitness  values 
and  assigned  probabilities  using  loumameni  selection.  Using  elitist  strategy  the 
suing  with  the  best  fitness  function  value  is  copied  to  the  new  population. 

3.  Produce  subsequent  generation  starting  from  the  cuircnt  using  crossover.  The 
strings  are  randomly  selected  for  mating  using  single  point  or  uniform  crossover. 
The  mating  between  two  strings  is  deterniined  by  the  probability  of  crossover  P,. 

A Perform  mutation  randomly  on  some  of  the  alleles  in  order  to  prevent  premature 
convergence  using  the  probability  of  mutation  P„. 

5.  Repeat  steps  3 and  4 until  convergence  is  achieved  or  the  maaimum  number  of 
generations  is  exceeded. 

The  simulation  optimization  procedure  is  externally  coupled  thus  first  a simulaiion 
is  run  and  the  results  that  form  the  decision  parameters  and  boundary  constraints  are 
extracted.  The  main  calling  program  has  a subroutine  that  reads  the  binary  and  text  files 
of  the  simulation  and  read  the  required  concentration,  head,  and  pumping  values.  In 
addition  to  the  values  read  from  the  simulaiion  result,  the  user  has  to  create  the 
parameters  required  to  run  the  core  operations  of  the  GA  module. 


TTie  objective  ftincKon  is  fotmulaled  using  the  ptu-amciers  and  constraints.  The 
optimization  proceeds  until  the  objective  function  is  satisfied  for  the  specified  population 
size  and  number  of  generations.  Once  the  optimization  results  are  obtained  the  new 
parameter  values  are  once  again  used  in  the  simulation  models  to  test  acceptability  of  the 
optimization  results  ui  terms  of  satisfying  the  management  objectives. 


GA  methodologies  in  current  use.  The  algorithm  was  designed  to  be  a simple  GA  with 
robust  working  mechanisms.  With  the  advent  of  newer  techniques,  individual  functions 
and  routines  ofihe  core  algorithm  inGENOPT  could  be  modified  in  the  future  to  create  a 
faster  and  more  efficient  search  algorithm. 

The  generation  ofrandom  numbers  is  the  first  process  in  the  algorithm.  In  order  to 
generate  random  numbers,  Knulh’s  algorithm  was  used  (Knuih,  1981)  as  described  in 
Press  etal.,  (1992).  The  random  number  generating  routine  returns  a uniform  deviate 
between  0.0  and  1 .0.  The  routine  works  by  reading  an  initial  seed  value,  after  that  a 
random  number  field  is  initialized  and  the  first  random  number  is  generated. 
Subsequently,  new  random  numbers  are  generated  using  a subhaclive  algorithm. 

Determination  of  the  population  size  was  based  on  the  average  size  of  schema  (the 
average  number  ofbitsperparameter(Goldberg,  1989).  In  a binary  coding  fonnulation 
of  strings  and  for  a schema  of  size  k and  string  of  length  I,  the  population  size  is 
determined  as  (Goldberg,  1992): 


where  ; is  Ihe  number  of  possibilities  for  each  chromosome  (cardinality  of 
chromosomes).  For  binary  formulation  the  value  of  is  2. 


The: 


I was  developed  on  the  basis  of  the  more  conventiotud 


(6-22) 


Encoding  of  the  decision  variable  imo  strings  (refeired  to  in  GA  as  chromosomes) 
was  accomplished  by  binaiy  coding.  It  is  also  possible  to  nse  integer  and  real  coding 
instead  of  binary  coding.  Form  decision  variables  in  an  optimisation  problem,  if  the 
decision  variables  are  binary  coded  as  binary  strings  of  length  n,  then  the  resulting 
chromosome  will  be  an  n x m string.  The  decoding  of  the  string  for  the  f decision 
variable  X;  was  performed  as  (Tang  and  Mays,  1998): 

+ (6-B) 

where,  a(i)  = (0,1 ) is  the  i**  binary  digit  value  in  the  corresponding  chromosome.  For 
example,  two  decision  variables  each  coded  as  binary  string  of  1 1 digits  will  have  a sling 
of  22  binary  digits  that  looks  like  lOOOIOI  1 101 101010001 10. 

The  fitness  of  a chromosome  was  evaluated  lo  determine  the  probability  that  a 
chromosome  will  be  selected  as  a parent  chromosome  lo  generate  new  chromosomes. 

This  was  achieved  through  the  use  of  a scaling  funclioo  5.  which  maps  objective  function 
values  to  an  appropriate  positive  fitness  value.  For  a vector  of  decision  variables  A"  the 
objective  function  value  g(A'(a»  Is  mapped  in  lo  the  fitness  value  /(a)  using  8; 

/(a)  = <5(g(X(a)))  (6-24) 

The  fitness  function  /(a)  that  determines  the  reproduction  of  a chromosome  is 
obtained  from  the  modified  objective  function  as  described  in  Equation  (6-21); 

g(A'(a)}=i(e,^,C.j..Ga„.) 

(6-25) 

For  the  fitness  function  described  in  Equation  (6-25),  the  scaling  function  was 
based  on  a linear  dynamic  scaling  (Tang  and  Mays,  1998).  This  method  ensures  non- 
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negaiive  fitness  function  for  all  chromosomes  and  the  worst  performing  chromosomes 
that  have  zero  fitness  will  be  excluded  from  reproduction.  The  linear  dynamic  scaling 
method  is  formulated  for  a given  chromosome  as: 

•f(/(a»  = £(e„.r„,eo„.)-Afm[/(a)(a,  sAtfr)]  (6-26) 

Selection  helps  choose  the  survivor  chromosomes  that  are  fit  for  the  new 
generation,  thus  forming  the  parent  chromosomes  for  the  next  generation.  For  the 
gencitiled  field  of  random  numbers  N.  thelolal  fitness  Miscalculated: 

S = S/(a.)  (6-27) 

where  /(fli)  is  fitness  of  the  i"  chromosome  in  the  old  population.  As  a next  step  the  k* 
parent  chromosome  is  selected  if 

s/(M  i/(o 

— <N£-^ — (6-28) 

A selection  scheme  was  implemented  based  on  the  fitness  of  the  individual 
chromosomes: 

(6-29) 

£/(»,) 

Crossover  makes  it  possible  for  alleles  of  chromosomes  to  undergo  changes  that 
will  allow  them  to  create  new  individuals  throu^  bit  combinations.  As  a result, 
crossover  is  considered  to  be  a very  important  search  operator  in  genetic  algorithms. 

Since  crossover  combines  useful  segments  of  different  parents  it  makes  it  possible  for  Ihe 
new  individual  to  benefit  from  advantageous  bit  combinations  of  both  parents  (Holland, 


1975).  Crossover  probQbiliiiesofO.6  lo  0.7  were  assigned,  and  depending  on  ihe  filness 
of  the  selected  parcnl  chromosome  crossover  was  performed. 

Mutation  helps  reintroduce  lost  alleles  into  the  population.  Alleles  are  feature 
values  assigned  to  chromosomes  in  natural  systems  or  strings  in  artificial  systems,  so 
each  string  is  assigned  an  allele.  The  reintroduciion  helps  bnng  back  bil  positions  lhat 
have  converged  lo  a certain  value  Ihroughoui  the  population  and  therefore  could  never  be 
regained  again  by  means  of  recombinalion  (Tang  and  Mays,  1998).  The  two  common 
forms  of  mutation,  creep  and  jump  mutation,  were  implemented  in  the  formulation  of  the 
solution.  Theoverallprobabiliiyofajumpmulalion  for  an  individual  Oj  and  number  of 
chromosomes  (binaiy  bits  in  a siring)  n,  was  deleimlned: 


Similarly  for  creep  mutation  with  number  of  parameters  the  creep  mutation 

probability  was  determined: 


For  overall  pipbability  of  mutation,  probabililies  of  jump  and  mulalion  and  creep 
were  set  lo  be  equal,  and  creep  mulalion  probability  was  expressed: 


If  Equation  (6-32)  is  expanded  using  binomial  expansion  and  ihe  lower  order  terms 
are  neglected  the  relationship  becomes  (Carroll,  1996): 


(6-30) 


(6-31) 


(6-32) 


(6-33) 


The  process  flowchart  of  ihe  GA  operations  is  presemed  in  Figure  6-2.  The 
externally  coupled  siniulation-opiiniiraiion  process  flowchart  is  presented  in  Figure  6-3- 


Figure  6-2.  Flowchart  of  core  operations  ofiheOA  algoritfun. 


Figure  6-3.. 


CHAPTER  7 

NUMERICAL  MODELING  RESULTS 


Applicaliod  ofCCGEN 

were  avsilabic.  The  global  sea  level  rise  scenarios  developed  by  the  IPCC  were  used 
■*lect  a mid-range  scenario  and  deiermioc  ihe  relative  sea  level  rise  for  the  specific 
lion  and  period.  The  results  are  presented  as  vulnerabilily  mops  and  a series  of  plots 


Research  conducted  in  different  parts  of  the  world  over  the  past  few  decades  clearly 
Mtes  that  coastal  areas  ate  under  the  increasing  threat  of  sea  level  rise.  In  the  United 
, the  state  of  Florida  is  one  of  the  states  with  a long  shoreline  (see  Figure  7-1 ) that 
the  possibility  of  being  impacted  by  future  sea  level  rise.  The  total  length  of 

and  banks  along  rivers  and  estuaries.  This  iong  coastline  makes  Florida  vulnerable 
impacts  of  climate  change  and  sea  level  rise,  especially  in  the  highly  populated 
tal  areas  that  face  an  increased  threat  due  to  the  combined  effects  of  both  climate 


ties  were  selected  tiom  south  Florida.  Digital  data  obtained  from  the  USOS  and 


t of  the  IPCC  to  generate  sea  level  rise  scenarios  for  Palm  Beach.  Broward,  and 
i-Dade  counties.  The  analysis  was  performed  using  the  mid  range  scenario  A I T 
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and  iha  respective  mean  sea  level  nends.  Maps  shovring  the  exleni  of  new  sea  level 
values  and  the  Lranslaticn  in  the  location  of  seacoasi  for  a simulation  period  of  choice 
were  generoied. 


A vulnerability  assessment  conducted  usings  digital  elevation  model  (DEM)  for 
the  entire  stale  of  Florida  made  it  possible  to  quantify  the  extent  of  possible  submergence 
due  to  sea  level  rise  (Tiruneb  and  Motz.  2001).  The  study  used  a 90  meter  cell  DEM.  Ihe 


proximity  to  sea  level.  The  result  showed  that  an  estimated  area  of  19,862  km^  is  at  sea 
level,  and  approximately  7.499  km*  has  an  elevation  of  1.0  m above  sea  level.  Table  8-1 


118 


DEM.  The  length  orcoe.<nline  for  each  of  these  counties  is  tabulaletl  in  Table  7-2. 


Tabic  7-2.  Length  of  coastline  for  counties  in  the  study  j 


studied  the  impnci  of  sea  level  rise  on  the  coastal  areas  of  Florida.  The  estimates  suggest 
that  sea  level  rise  along  the  Florida's  shoreline  has  accelerated  at  ante  of  approximately 
20  to  40  centimeters  (g  to  1 6 inches)  per  century  since  1 932.  This  figum  is  more  than  six 
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times  the  rale  recorded  by  earlier  tide-gage  records  and  estimates  from  the  geological 
history  of  the  past  three  ihousand  years.  The  relative  sea  level  rise  along  the  shorelines  of 
Florida  is  presented  in  Table  7-3.  Relative  sea  level  accounts  for  fiiciors  such  as 


shoreline  degradation,  subsidence,  and  other  geological  inslabilities.  and  is  therefore 
larger  than  euslalic  sea  level. 


Table  7-3.  Estimated  rates  of  relative  sea  level  rise  for  selected  sites  in 


Femandma  Beach 


The  Center  for  Operational  Oceanographic  Products  and  Services  (CO-OPS)  is  the 
part  ofthe  National  Oceanic  and  Atmospheric  Administration  (NOAA)  responsible  for 
collecting,  analyzing  and  distributing  historical  and  real-titne  observations  and 
predictions  of  water  levels,  coastal  currents  and  other  meteorological  and  oceanographic 
data.  The  water  level  records  registered  by  the  CO-OPS  stmions  are  a combination  of  the 
fluctuations  ofthe  ocean  and  the  vertical  land  motion  at  the  location  ofthe  station.  The 
sea  level  variations  determined  are  the  linear  trend,  average  seasonal  cycle,  and  inter- 
annual variabili^  at  each  station.  Mean  sea  level  trends  for  selected  stations  in  south 
Florida  are  presented  in  Table  7-4.  Monthly  data  up  to  the  end  of  1 999  were  used  in  the 
calculation,  and  all  slalions  had  data  spanning  a period  of2S  years  or  more. 
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since  the  1930s.  When  the  observed  local  tread  is  combined  with  the  euslnlic  sea  level 
rise  estimated  by  Global  Circulation  Models,  the  corresponding  estimates  ofpmbabic 
relative  sea  level  rise  ore  obtained.  The  mean  monthly  sea  level  observations  collected 
by  the  CO-OPS  stations  in  Florida  ore  presented  in  Figure  7-5.  The  values  were  averaged 
monthly  from  houriy  tide  records  observed  at  the  respective  stations.  The  mean  sea  water 
level  is  defined  with  respect  to  the  Mean  Lower  Low  Water  (MLLW).  The  MLLW  is  a 
tidal  datum  defined  by  the  average  of  the  lower  low  water  height  of  each  tidal  day 
observed  overthe National  Tidal  Datum  Epoch.  For  stations  with  shorter  series, 
simultaneous  observational  comparisons  are  made  with  a control  tide  station  in  order  to 
derive  the  equivalent  datum  ofthe  National  Tidal  Datum  Epoch.  The  (NTDE)  is  a 
specific  1 9-year  period  over  which  tide  observations  are  taken  to  determine  Mean  Sea 
Level  and  other  tidal  datums  such  as  Mean  Lower  Low  Water  and  Mean  High  Water.  The 
latest  update  defines  the  1 9-year  period  as  1983-2001.  A tidal  datum  is  a veitir»l 
reference  based  on  a specific  stage  of  tide,  which  serves  as  a bn.scline  elevation  to  which 
sounding  depths  or  topographic  heights  are  referenced.  The  1 9-year  period  includes  an 
1 8.6  year  astronomical  cycle  that  accounts  for  all  significant  variations  in  the  moon  and 
sun  that  cause  slowly  vajying  changes  in  the  range  of  tide.  The  NTDE  has  been  adopted 
so  that  tidal  datum  determinations  throughout  the  United  Slates  will  be  based  on  one 
specific  common  reference  period. 
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Figurr  7'S.  Monthly  Mean  Sea  Level  plot  for  stationa  in  South  Florida. 

Ihe  monthly  mean  sea  level  plot  presented  in  Figure  7-5  shows  that  there  is  a 
gradual  rise  for  all  the  stations  through  the  years  of  record.  Comparatively,  the  highest 

level  rise  for  select  stations  were  performed  using  the  trend  values  and  all  scenarios.  The 
results  presented  in  Figure  7-6  show  a gradual  rise  of  about  0.4  to  0.6  m depending  on  the 
specific  scenario  selected.  A middle  of  the  range  value  is  obtained  by  scenario  AIT. 
which  has  a reasonable  probability  of  occurrence 

The  study  conducted  and  the  results  oblairted  clearly  demonstrale  that  sea  level  rise 
will  result  in  enhanced  coastal  erosion,  coastal  flooding,  loss  of  coastal  wetlands,  and 
increased  risk  from  storm  surges,  particularly  in  Florida  and  much  of  the  U.S.  Atlantic 


the  uitique  topography  of  the  coastal  areas  of  Florida  and  the  Allanlic  coast. 
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it  ha5 10  be  noted  (hat  the  development  oflhe  IPCC  scenarios,  which  is  based  on  the 
SRES,  does  not  places  likelihood  or  probability  of  occurrence  on  any  of  Ihe  scenarios.  It 
is  therefore  suggested  that  any  research  to  investigate  the  impacts  of  sea  level  rise  should 
recognize  the  inheient  uncertainties  of  long-term  projections.  While  these  projections 
provide  a long-temi  context  for  Ihe  analysis  of  a change  expected  to  happen  in  the 
relatively  near-term,  the  deuiiled  analysis  should  consider  the  localization  of  the  problem 
in  greater  details.  Asa  result,  the  projections  presented  in  the  previous  plots  and  in  Table 
7-S  should  not  be  inteipreled  as  ~whal  is  bound  to  happen."  Instead,  they  should  be 
inicrprcted  as  "what  could  happen."  Moreover,  the  calculation  of  shoreline  tninslalion 
has  to  be  done  at  a much  smaller  scale  in  order  to  account  for  the  existence  of  localized 
infra.ciructures  and  mitigation  impacts.  The  shoreline  translation  presented  in  Table  7-S 
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2025  I 2050  I 2075  | 2100  | 2125  | 2150  | 2175  | 


).1S2  I 0.281  I 0.423  | 0.576  | 0.742  | 0.919X1 


Applicalian  ofSEAWATSLR 

The  modified  seawmer  intrusion  code  was  used  to  model  groundwater  How  in 

11ie  models  include  a cross<scc(ional  model  from  south  Florida  and  a three'dimensional 
model  from  the  Silifke  Della  in  Turkey. 

Due  to  the  volume  of  work  involved  primarily  doe  to  the  transient  nature  of  the 

was  used  in  the  simulations.  The  AIT  scenario  that  represents  the  mid-range  was  used 
for  all  the  cases.  The  results  presented  in  this  chapter  show  a comparison  of  simulations 
attheendofTOyearsand  100  years  from  present.  The  transient  simulations  in  between 
presented  in  the  appendix. 


Since  analysis  of  this  nature  primarily  focuses  on  changes  in  pans  of  the  aquifer 
subject  10  saltmuer  intrusion,  it  becomes  necessary  lodeijne  a parameter  that  describes 
llie  spatial  extent  of  changes  in  concentration  for  Imnsieni  sitnulalions.  For  the  purpose 
of  this  analysis,  a parameter  here  after  referred  to  as  zone  of  intluence  is  used.  A zone  of 
influence  is  defined  as  part  ofthe  aquifer  that  will  be  subject  to  significant  changes  in 
terns  of  aquifer  salinity  due  to  transient  changes  of  stresses.  The  zone  of  influence  is 
characterizedby  a dimensionless  aquifer  salinity  change  ratio.  The  difTcrence  in 

monitored  concentration  at  a given  location  x between  rime  steps  i and  f+l  is  given  as: 


while  the  dimensionless  number  ofthe  zone  of  influence  is  given  as: 


(7-2) 


Cross-Sectional  Model 

The  cross-sectional  model  was  originally  developed  to  simulate  the  groundwater 
discharge  patterns  near  the  coasiof  Biscayne  Bay  (Langevin,  2001).  The  study  area  is 
pan  ofthe  surficial  aquifer  system  in  south  Florida.  A briefdescription  ofthe 
hydrogeology  and  stratigraphy  documented  m the  report  (Langevin,  2001)  is  presented 
here  (ogive  background  information  ofthe  area. 

South  Florida  presents  unique  hydrologic  feaiures  primarily  due  to  the  dynamic 
inietaciion  between  surface  water  and  groundwater.  Rainfall  and  evapolranspiration  in 
the  area  are  generally  hi^  with  marked  seasonal  variability  of  rainfall.  For  the  study 
area,  the  average  rainfall  for  the  period  from  Januaiy  1 989  to  September  1 998  was  eboul 
141  centimeters  per  year.  About  75  percent  ofthe  rainfall  was  recorded  during  the  wet 
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season  that  enends  from  June  lo  October.  Direct  mca.surements  of  evapolmispiration  in 
parts  of  South  Florida  resulted  in  a range  of  122  to  130  centimeters  per  year.  AJthough 
the  measurements  may  not  be  representative  of  the  entire  southern  region  of  Florida,  the 
talcs  of  ET  are  as  high  as  90  percent  of  the  rainfall  in  some  ports. 

The  hydmstiatigrapby  of  souUiem  Florida  is  primarily  composed  of  the  shallow 
surficial  aquifer  systems  and  the  deeper  Floridan  aquifer  system.  As  observed  from  the 
water  table  maps  presented  in  the  report  (Langevin,  2001),  the  general  flow  of  ground 
water  is  toward  the  coast,  thus  implying  that  the  primary  hydrologic  component  that 
enables  the  aquifer  lo  maintain  the  water  table  is  net  recharge  from  rainfall.  The 
Biscayne  aquifer,  which  is  pan  of  the  surficial  aquifer  system,  consists  of  parous 
limestone  from  the  Pliocene  to  Pleistocene  epochs.  The  stratigraphic  map  of  Florida 
created  from  data  obtained  from  Florida  Geological  Survey  is  presented  In  Figure  7<I7. 
A map  showing  the  relations  of  geologic  formalians,  aquifers,  aud  semi  permeate  units 
of  the  surflcial  aquifer  system  across  ccnlml  Miami-Dade  County  is  presented  in  Figure 
7-18  (Source  SOFIA  website  hlln:/,'siUia.use.s.covinuhlicalions/wri,'0ll-l75n 
Coconut  Grove 

The  Coconut  Grove  model  modified  as  pan  of  this  research  represents  iJie  coastal 
aquifer  south  of  Miami  in  Miami-Dade  County.  The  cross-section  extends  for  5.900  m 
Inland  and  2.100  m towards  die  sea.  The  base  of  the  model  extends  lo  40  m below  sea 
level.  For  the  model  discretmilion,  a uniform  column  width  of 200  m and  uniform  layer 
thickness  of2m  were  used.  The  model  has  a constant  head  boundary  on  the  seaside  and 
a constant  flux  boundary  on  the  inland  side.  The  freshwater  head  for  the  constant  head 
boundary  cells  was  calculated  from  a downstream  stage  value  of  0.22  m,  which 
tvpn!scm.s  an  average  value  for  sea  level  from  1989  lo  1998.  A net  recharge  value  of  38 
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cm/yr  was  assigned  la  the  overlying  recharge  boundary.  A conslani  seawater 
concentration  ofJSkg/m’ was  used  for  the  seaside  boundary.  The  model  calibration  data 
was  obtained  from  Ihe  USGS,  and  on  the  basis  of  the  calibration  results,  the  following 
parameters  were  used  for  the  simulation. 


of  Florida  and  Florida  Keys. 
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Table  7-6.  Model  pattmawt  mad  for  Coconul  Grove  simuUlion 


Mode]  Parameaer 

Unit 

Value 

K.  niorizcciial  hvdriulie  conductfvltvt 

rtitXav 

9000 

1 

1 

1 

1 

,8 

9 

(b  tinneiojdinal  rlisnaesjvirvt 

1 

s 

1 

1 

1 

0 1 

n tnomsltvt 

dlmenmosikaa 

OJ 

(Vm  Ifliixl 

m'/d/m 

15 

R fner  rerharoet 

m/VT 

0 3g 

dimenskmless 

0.2 

Figure  7‘I9.  Discretieation  of Cocomji  Grove  model. 

In  orderioesiablish  a more  realiaiic  inillal  condillon,  a simulation  was  run  for  one 
million  years  using  SEAW AT.  Performing  the  simulation  for  such  a long  period  made  it 
possible  for  the  system  to  ansin  a stale  of  dynamic  equilibrium.  The  hydraulic  head  and 
concentration  values  at  the  end  of  the  million  years  simulnlion  were  used  as  inilial 
conditions  to  perform  simulations  for  the  next  1 00  years.  Simulations  were  tun  forarolal 
of  100  years  assuming  uniform  recharge  and  flux  values  to  persist  throughout.  The 


simulation  using  SEAWATSLR.  The  AIT  sea  level  rise  scenario  was  assumed,  which 
resulted  in  a sea  level  rise  of  about  OJ  m by  the  year  2100.  The  transport  compcmenl  was 
solved  using  the  method  of  characteristics.  The  results  show  a change  in  the  position  of 
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ihc  iransilion  zone  due  toihe  incluaicn  of  a sea  level  rise  componeru  in  the  consunr  head 
boundaty.  Monitoring  wells  weie  introduced  along  the  transition  zone,  and  hydrogtaphs 
of  concentration  were  plotted  to  see  the  iiansieni  change  in  the  arnbient  conceniralion. 
The  simulations  were  performed  with  steady  values  of  parameters  and  no  additionnl 
pumping  stresses.  The  purpose  of  this  nppronch  was  to  see  the  le^onseofthe  system  to 
natural  conditions. 


Figure  7<20.  Location  of  monitoring  wells. 

The  monilored  concentration  values  plotted  in  Figure  7-21  eahibii  changes  in 
concentration  due  to  sea  level  rise.  The  changes  as  presented  in  Figure  7-22  vary  with 
the  location  of  the  monitoring  wells.  The  general  trend  is  that  the  rate  of  change  tends  to 
decrease  as  the  monitoring  location  is  farther  inland. 

In  order  to  have  a more  realistic  inhiai  condition,  a no-sea  level  rise  simulation  of 
5.000  yeais  performed.  The  final  concentration  and  head  values  were  used  to  define  the 
initial  condition  for  the  1 00  yean  simulation  of  the  period  until  the  year  2100.  The 
^liles  of  concentration  were  plotted  for  the  entire  period  of  simulation  at  specific  time 
intervals  by  exporting  the  data  to  Surfer^  and  krigging  to  create  contours.  In  order  to 
identify  Ihc  areas  that  will  be  subject  to  significant  changes  due  to  variability  of  stresses, 
the  zones  of  influence  were  also  ploued  for  the  same  time  periods.  The  results  show  that 
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thu  portion  of  ihc  aquifer  that  will  be  subject  to  significant  changes  due  to  sea  level  rise 
under  natural  conditions  is  the  area  that  extends  from  3.000  m to  about  4,200  m inland. 
The  areas  extending  front  the  shore  up  to  3,000  m inland  will  be  subject  to  moderate 
changes  in  concentration  under  the  conditions.  The  practical  implication  of  this  anaijsis 
is  that,  the  better  location  of  pumping  wells  from  the  seawater  intrusion  point  of  view 


Figure  7-21.  Plot  of  monitored  concentrations  for  Coconut  Grove  monitoring  wells. 


Figure  7-22.  Plot  of  percent  changes  in  monitored  concentration  due  to  .sea  level  rise. 


•TO  3m  4000  j3oo  ^i5o  ! 

CiBacceOnO 

Figure  7'23.  Profile  arccnceotralion  after  100  years  for  Coconut  Grove. 


Figure  7-24,  Zone  of  Influence  aPerl  00  years. 

In  ihe  plots  of  concentration  From  the  monitoring  wells,  there  Isa  noticeable 
fluctuation  in  values  over  the  simulation  period.  In  order  to  understand  Ihe  phenomena, 
groundwater  flow  velocity  and  discharge  were  plotted  for  Ihe  cross-section.  The  results 
show  that  groundwater-seawater  flux  exchange  occurs  underneath  the  coast  resulting  in 


Iheobserved  oscillations  in  monitored  concentration  from  wells  close  to  Ihe  coast.  The 


observed  osciliailons  tend  to  dampen  as  the  wells  are  located  farther  away  from  the  coast. 


Accordingly,  well  1.  which  is  the  closest  to  the  surface,  exhibits  Ihe  highest  oscillation. 


while  well  4 exhibits  the  least  oscillation.  The  groundwater  discharge  generally  r 
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The  comours  of  head  plonedin  figure  7-27  showen  increase  duelo  the  sustained 
SCO  level  rise.  The  gradient,  however,  still  follows  the  same  trend  with  decreasing  head 
values  towards  the  coast.  In  the  scenario  where  sea  level  rise  was  not  considered,  the 
contours  tend  to  rise  with  a mild  slope  after  the  active  zone  of  exchange,  while  in  the  sea 
level  rise  scenario,  the  rise  appears  to  be  sli^tly  higher.  This  increase  shows  the 
response  of  the  constant  head  boundary  and  its  significance  in  Influencing  the  flow 


Figure  7-27.  Profile  of  head  at  the  end  of  IDO  years  simulation. 

Three-Dfisieiislonal  Model 

The  three-dimensional  model  presented  in  this  study  represents  s deltaic  aquifer 
close  10  the  Mediterranean  coast  of  Turiey.  The  available  data  for  this  model  was  very 


aquifer.  The  primary  form  of  aquifer  data  was  generated  from  digitized  topographic 
maps  and  hand  sketches  of  cross-sections  obtained  from  geologists  who  have  some 


knowledge  of  the  area.  Some  data  on  water  use,  aquifer  water  conductivity,  and  river 
stage  were  also  obtained.  Although  the  data  were  sparse  and  insufficient,  care  was  taken 
to  represent  the  current  condition. 

The  conductivity  data  were  convened  to  total  dissolved  solids  (TDS)  using 
appropriate  regression  equations.  A two  dunensional  simulation  of  the  Coksu  Della 
aquifer  was  performed  earlier  (Gordu,  2000).  The  study  included  deleiminaiion  of  some 
aquifer  patamclcrs  that  helped  in  the  development  of  the  three  dimensional  model  used  In 
this  study. 

Goksu  Delta  at  SIlifke 

The  Goksu  Della  is  predominantly  an  alluvial  deposit  formed  by  the  Goksu  River, 
which  is  about  260  km  long  and  has  a drainage  area  of  10,400  km'.  The  Delta  has  an 
area  of  about  ISO  km’  with  a semi-arid  Mediterranean  climate  (Gordu,  2000).  The  Goksu 
River  with  an  annual  average  discharge  of  about  9.5  million  m’/day  is  a significant 
component  of  the  hydrology  of  the  della.  There  are  two  lakes  of  significance  in  the  della, 
namely  Akgol  and  Pomdenia.  Akgolisa  freshwater  lake  with  on  area  of  oboui  12.45 
million  m’  while  Paradeniz  is  a saltwater  lake  with  an  area  of  about  4.05  million  m’. 

There  is  0 man-made  canal  that  connects  the  two  lakes.  The  area  is  an  agriculniral  area 
with  waler  demands  for  bolh  irrigation  and  water  supply  of  the  population. 

The  primary  source  of  geologic  data  used  by  Gordu  (2000)  and  the  current  case 
study  is  the  hydrogeologic  iovostigalion  conducted  by  Haar  and  Heunks  (1992).  The 
results  of  their  investigation  included  lithological  piofiles  from  which  it  was  concluded 
that  the  coarse  grained  medium  ranging  from  50  to  1 30  m in  depth  is  the  aquifer  with 
connections  to  the  adjacent  karstic  formations  and  inlemediaie  ihui  layers  of  clay. 
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The  invesliiadon  siso  included  coUecling  bodi  surface  and  ground  woler  samples 
from  which  a woler  table  map  was  constnicled  fi3r  the  surficial  phreatic  aquifer.  Based 
on  the  map,  there  is  a hydraulic  gradient  towards  the  sea  where  the  water  table  is  almost 
at  sea  level,  nie  deeper  aquifers  are  confined  by  an  overlying  scmi-permeoble  clayey 
confining  zone.  It  is  also  mentioned  by  Gordu  (2000)  that  Mencngic  (IW!)  conducted  a 
study  about  the  potential  risks  of  groundwater  pollution  in  the  Goksu  della.  On  the  basis 
of  the  results  documented  by  Haar  and  Heunks  (1992)  and  Menengic  (1998).  Gordu 
(2000)  presented  a stratigraphic  classification  of  the  Goksu  Delta.  The  classification  is 


adapted  in  this  study  and  presented  in  Table  7-7. 
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Previous  studies  have  estimated  the  hydraulic  conductivities  to  range  liom  4 to  500 
in/day  in  the  horizontal  direction  and  Cnjin  0.05  to  27  m/day  in  the  vertical  direction,  and, 
thus,  the  aquifer  is  relatively  homogeneous  and  anisotropic.  The  dispetsivity  of  an 
aquifer  is  a parameter  known  to  be  dependent  on  the  scale  of  the  flow  model  length  along 
which  the  parameter  is  defined.  Depending  on  the  geometry  of  the  cross-section,  values 
of  1 000  m and  50  m were  used  respectively  in  the  study  hy  Qordu  (2000). 

The  spatial  discretizalion  of  the  model  is  based  on  a finite  difference  grid  of  42 
rows,  3g  columns  and  20  layers.  The  spacing  of  the  columns  is  kept  at  500  m while  the 
layer  thickness  is  variable  with  maximum  depth  extending  to  700  m below  sea  level.  The 
boundaries  include  a constant  head  boundary  with  constant  sea  water  eonccniiation  of 
39.4  kg/m^  (representative  of  the  Mediterranean  sea),  a no- flow  boundary  that  extends 
from  northwest  to  southeast  of  the  model  area  defining  the  bedrock,  a general  head 
boundary  for  the  inland  portion  of  the  aquifer  in  the  northwest,  lakes  and  a network  of 
drains.  The  stage  for  the  river  and  drains  was  set  from  elevation  and  stage  data  available. 
For  the  purpose  of  detailed  analysis,  monitoring  wells  were  strategically  located  close  to 
the  boundary  and  surface  water  features.  The  constant  head  boundary  Is  set  to  extend  for 
500  m everywhere.  Perhaps  with  the  availability  of  bathymetric  data  in  the  future,  this 
boundary  will  have  to  be  specified  differently.  At  this  stage,  the  model  is  based  on  a 
coarse  grid,  which  may  have  to  be  modified  in  the  future  with  the  availability  of  more 
geophysical  and  hydrologic  data.  The  spatial  discretization  of  the  model  is  presented  in 
Figure  7-28.  and  the  model  boundary  features  are  presented  in  Figure  7-29.  The  model 
parameters  used  in  the  Silifke  model  are  presented  in  Table  74- 


In  order  w analyze  ihe  resulls  of  aimulaiion  for  the  Goksu  Delia  aquifer,  two  cross- 
seciions  were  considered.  The  first  cross-secUon  (seclion  Y-Y)  was  along  row  20  with  an 
easi-west  orienlollon,  while  ihe  second  cross-seclion  was  along  column  1 9 (section  X-X) 


Unit 

Value 

(hnrironul  hydraulic  conduclivitvl 

m(dav 

200 

K.  tvenical  hvdnuiic  conduciivilv) 

0.05 

100 

100 

dimensionleas 

0.1 

0 

m/vr 

0-1.00375 

S.  (specific  yield) 

dimensionless 

01 
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The  simulaiion  results  are  presenced  in  the  fpnn  afcrosS'Seclipnal  plots  along  both 
directions.  Zones  ofinflijence  were  also  constructed  from  the  crosS'Scctions  considered. 
Plots  of  monitored  concentration  versus  time  are  also  presented  to  see  the  evolution  of 
concentration  through  time  at  a given  location. 

Contours  of  equivalent  freshwater  head  are  also  plotted  for  the  cross-sections 
considered  and  the  transient  variability  is  compared  for  both  scenarios.  Contours  of 
equivalent  freshwater  head  are  also  plotted  for  the  ^rst  layer.  This  makes  it  possible  to 
observe  the  gradient  of  the  groundwater  head  and  make  some  observations  regarding  the 


flow  direction. 


Figure  7-30.  Conccntraiion  profile  comparison  along  X-X  after  100  years. 
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Figure  7-3 1 . Concenirelion  profile  comparison  along  Y-Y  after  100  years. 


iCoksu 
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ccndilions,  Ihe  lower  aquirens  tend  to  reach  a limiting  value  of  concentration  explains 
partially  why  the  difference  in  concentration  decreases  vertically  downwards.  There  are 


deposition,  but  ntimerical  dispersion  seems  Ihe  more  likely  cause.  Similar  phenomena  of 
localized  increase  were  observed  in  the  simulation  results  of  the  Burdekin  Della  in 
Australia*.  The  simulations  tor  the  Burdekin  Della  used  ihe  two'dimensionol  version  of 
SUTRA. 

One  interesling  observation  from  Ihc  simulated  concentration  profiles  is  die 
significance  of  the  geomeliy  of  the  cross-section  in  influencing  the  shape  of  the  transition 
zone  profile.  Due  to  die  existence  of  a no  flow  boundary  with  gradually  decreasing 
extent,  the  shape  of  the  transition  zone  exhibits  a curvature  that  follows  Ihe  geomeuy  of 
the  cross-section. 

As  a result  of  Ihe  combination  of  the  hydrodynamic  and  variable  density  pioccsses 
induced  by  sea  level  rise,  parts  of  the  aquifer  that  are  expected  to  undergo  significant 
changes  in  salinity  are  presented  in  Figure  7-35.  Considering  section  X-X,  the  zone  of 
influence  extends  from  approximately  4.000  m from  Ihe  seaside  boundary.  With  an 
average  thickness  of  about  4,000  to  5.000  m.  it  extends  to  about  400  m depth  and  starts  a 
reversal  due  to  Ihe  geometry  of  the  no  [low  zone.  Considering  section  Y-Y.  the  zone  of 
influence  extends  from  about  1,500  Inland  to  about  1 1,000  m where  it  expands 
significantly  in  the  middle  portion  of  Ihc  aquifer  (about  200  to  400  m depth)  and 
continues  to  thin  out  towards  the  west.  Therefore,  the  influence  of  sea  level  rise  on  a 
century  scale  for  Coksu  Della  is  not  expected  to  be  sigoificant  fanhei  than  about  5,000  m 
along  X-X  and  10,000  m along  Y-Y. 


: the  maximum  specified  in  Ihe  CHD  suggesting 


. KumarNsTaysfiafCStRO,  AustTnlia  November  2002 
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From  the  contour  of  groiindwotcc  observed  from  each  layer,  the  general  flow  trend 
Is  towards  the  sea  with  appreciable  gradient  oriented  from  northwest  to  southeast.  This  is 
an  indicative  of  groundwater  discharge  towards  the  Mediienanean  Sea.  In  order  to 
confimt  this  theory  it  would  be  necessary  to  investigate  the  seasonality  of  the  contours 
usbtg  data  obtained  from  a network  of  monitoring  wells.  Also,  from  the  groundwater 
contour  of  layer  1,  the  Goksu  River  appears  to  be  a generally  gaining  river.  The  results 
presented  are  to  be  used  as  preliminary  simulation  results  obtained  from  a sparse  dataset, 
and  therefore  (he  Goksu  Della  model  needs  to  be  refined  in  the  fliture  with  more 
coljbralion  work  os  additional  data  sets  become  available. 

Simulated  Heads  and  Conccotratious 

The  changes  in  simulated  head  and  concentration  values  of  a coastal  aquifer  system 
are  the  primary  forms  of  the  response  to  a rismg  sea  level  boundary.  The  precise  nature 
of  the  response  could  be  described  by  a proper  characierizalion  of  the  seaside  boundary. 
General  observations  from  the  simulations  perfomted  make  it  possible  to  ctsncludc  that 
there  is  appreciable  increasing  trend  in  both  simulated  head  and  concentration  due  to  sea 
level  rise. 

In  the  current  formulation  of  the  constant  head  boundary  component,  the  monitored 
head  values  are  noted  to  rospond  to  the  linear  increase  of  sea  level  rise.  However,  if  the 
inctease  m sea  level  is  formulated  to  have  an  exponential  trertd,  the  response  is  likely  to 
be  slightly  dilTerent-  The  difference  in  head  values  between  a no  sea  level  rise  scenario 
and  the  one  that  accounts  for  this  is  a factor  of  the  rate  of  sea  level  rise  incorporated  in 
the  simulation.  Since  the  sea  level  rise  scenario  simulations  are  based  on  a futuristic 
assumption  based  on  a number  of  factors,  it  is  not  possible  to  calibrate  the  model  for  a 
given  scenario.  It  could  also  be  argued  that  the  linear  increase  in  monitored  head  values 


could  have  a difTerent  characlorisdc  by  uufoducing  a distance  dependent  propagation 
factor.  One  way  of  detennining  this  factor  would  be  to  conduct  a controlled  laboratory 
experiment  where  a static  saltwater  water  column  serves  as  a seaside  boundary.  This 
setup  issimilario  the  analog  models  used  by  researchers  a few  decades  ago.  Tbe  results 
of  this  type  of  experiment  will  help  determine  the  variation  in  monitored  heads  in 
response  to  arise  in  the  sialic  water  head.  By  applying  appropriate  similitude  analyses, 
the  results  could  he  used  to  further  refine  the  numerical  model  developed  in  this  study.  In 
any  case,  the  fact  remains  that  the  lime  variant  constant  head  boundary  plays  a signiftcani 
role  in  influencing  the  flow  regime  and  simulation  values. 

Shift  of  the  Seawater-Freshwater  Interface 

Under  steady  natural  stress  conditions  where  there  is  no  fteshwaler  withdrawal,  the 
shifi  of  the  seawater-freshwater  interface  is  brought  about  by  the  dynamic  interaction 
between  the  seawater  and  freshwater.  The  primary  driving  force  in  the  vicini^  of  the 
coast  appears  10  be  the  advective  flux  as  it  is  observed  in  the  flow  velocity  and 
groundwater  discharge  vectors.  Traversing  farther  inland  along  a cross-section,  the 
dispersive  component  begins  to  lake  effect  lesulling  in  appreciable  spread  in  the 

The  impacts  of  Sea  Level  Rise  on  Groundwater  Flow  Regimo 

The  general  impact  of  sea  level  rise  on  the  groundwater  flow  regime  is 
chaiacteriaed  by  the  creation  of  a very  active  mixing  zone  dose  lo  the  coast.  In  this  zone, 
seawater  enters  at  the  bottom  of  tbe  aquifer  driven  by  the  hydrostatic  piessure  variation 
and  the  density  difference.  After  mixing  with  the  ambient  aquifer  water,  the  general 
hydraulic  gradient  results  in  discharge  of  water  towards  the  sea.  In  the  long  run,  this 
process  tends  to  shift  the  inier&ce  further  until  a state  of  dynamic  equilibrium  is  reached. 
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Since  squIfcrsuvaUo  subject  loadditicnBl  stresses  such  as  puinplng.  it  is  rather  difficult 
tosa>  exactly  how  long  it  takes  the  aquifer  to  reach  this  state  of  dynamic  equilibrium. 
Results  obtained  Irom  a simulation  performed  fora  period  of 3.fX)0 years  were  plotted 
and  analyzed.  The  plots  ofmonitoned  concentration  versus  time  exhibited  a steady  slate 
after  approximately  2,500  years  Ibr  most  of  the  layers.  Even  after  3,000  years,  there  are 
portions  of  the  aquifer  that  have  not  attained  dynamic  equilibrium  as  it  is  seen  in  layers 
6,7,0,  and  1 1.  Ba.sed  on  this  observation,  it  appears  that  under  natural  constant  stress 
conditions,  the  aquifer  Is  likely  to  take  much  longer  than  3,000  years  to  reach  the  stale  of 
dynamic  equilibrium.  When  the  impacts  of  a dynamic  sea  level  change  and  pumping  are 


included,  the  variation  in  aquifer  stress  is  aggravated,  thus  increasing  the  time  required  to 


attain  dynamic  equilibrium. 


Figure  7*37.  Monitored  concentration  at  well  10  after  3000  years  simulation. 

The  sensitivity  ofaquifers  to  sea  level  rise  and  the  impact  thereof  ore  governed  by  a 
number  of  factors.  One  such  factor  is  geometry  of  the  aquifer.  Aquifers  with  smaller 
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geometries  are  likely  to  uodergo  dynamic  changes  relatively  faster  than  aquifers  of  larger 
geometry.  The  second  factor  Is  the  inherent  aquifer  properties  primarily  hydraulic 
conductivity  and  transmissivity.  Aquifers  with  smaller  transmissivity  values  are  likely  to 
take  longer  to  respond  to  changes  in  sea  level. 

The  time  required  for  a semi'confined  or  confined  aquifer  as  compared  to  an 
unconfined  aquiferis  relatively  shorter  in  small-scaie  modeling.  In  a regional  scale 
modeling,  this  response  lime  is  generally  within  the  scale  of  centuries.  In  effect,  ihe 
impact  of  a sea  level  rise  ofO.5  m is  fell  after  many  decades.  Even  the  impossible 
condition  of  a sudden  sea  level  rise  is  not  likely  to  result  in  sudden  rise  in  aquifer  salinity 
due  to  the  nature  of  groundwater  flow  velocity.  The  actual  groundwater  flow  velocity  is 
governed  by  the  porosis  of  the  medium  and  porosity  also  plays  a significant  role  in  the 
process.  In  coastal  aquifers  where  there  is  a possibility  of  secondary  porosity  due  to 
fracturing  or  if  there  are  sedimentary  deposits  where  localized  reactions  result  in  solution 
channels  where  the  flow  may  not  necessarily  be  governed  by  Darcy’s  law.  the  process  is 
expected  to  be  accelerated. 

Based  on  the  observations  from  the  various  simulations  performed  to  underhand 
Ihe  process  of  saltwater  intrusion  and  the  factois  that  govern  it,  generalized  quaillaiive 
observations  were  made.  The  observations  are  plotted  in  Figure  7<3S  to  describe  the 
significance  with  lime.  The  impacts  arc  classified  as  high,  moderate,  ur  low  and  Ihe  time 
scale  is  presented  in  a scale  ranging  from  century  to  millennia.  The  figure  depicts  the 
trend  of  response  in  terms  of  aquifer  salinity. 


Century 


Mill  enium 


KiUeiua 


Figure  7>3g.  Schemotic  presenuiioit  of  (rend  of  response  of  coastal  aquifers  to  stress 
The  interpretation  of  the  qualitative  description  presented  hi  Figure  7<38  is  as 


a century  scale.  However  if  pumping  is  continued  at  a given  rale  for  a longer  time  the 
aquifer  will  eventually  reach  the  slate  of  dynamic  equilibrium  on  a lime  scale  of 
millennia  at  which  lime  the  rate  change  becomes  smaller.  Even  If  pumping  is  increased 
Ihroughouu  Ihere  will  be  a point  where  the  aquifer  becomes  completely  saline  or  dry  as 
the  case  may  be  and  the  pumping  will  have  no  addillonal  impact  to  measure. 

Aquifers  with  higher  transmissivity  and  hydraulic  conductivity  values  are  likely  to 
attain  a stale  of  dynamic  equilibrium  within  centuries.  On  the  other  hand  aquifers  with 
smaller  Iransmissiviiy  values  might  take  millennia  to  do  the  same.  Porosity  of  an  aquifer 


dynamic  equilibrium.  Bui  on  a scale  of  millennia  it  is  likely  to  have  a moderate  elTeci  in 


159 

inducing  dynamic  equilibrium  especially  if  a Lime  scale  is  long  enough  for  facTors  such  as 
secondary  porosity  lo  lake  effect. 

The  geomelry  of  the  aquifer  has  a low  bur  significant  influence  on  the  state  of 
dynamic  equilibrium  on  a time  scale  of  a century.  As  lime  scale  increases,  this  influence 
also  increases.  For  example  if  two  aquifers  of  the  same  hydraulic  parameters  are  set  lo 
have  dificrcnl  geometries  where  one  is  larger  and  the  other  is  smaller,  then  the  smaller 
geometry  aquifer  is  likely  to  respond  fo-slcr  in  terms  of  dynamic  equilibrium  compared  to 
the  larger  geometry  aquifer. 

Model  Sensitivity  Analysis 

Sensitivity  analysis  performed  to  evaluate  the  signilicancc  of  the  various  aquifer 
parameters  and  boundary  conditions  makes  it  possible  to  see  how  the  various  parameters 
influence  the  flow  regime.  The  parameters  used  in  the  sensitivity  analysis  were  hydraulic 
conductivity,  flux  boundary,  dispersivity,  recharge,  and  pumping.  The  cross'sectional 
models  developed  for  the  sensitivity  analysis  were  all  variations  of  the  surficial  aquifer 
model  explained  earlier.  For  the  pumping  sensitivity  enalysis.  five  pumping  wells  were 
located  at  equal  distances  from  each  other  as  shown  in  Figure  7-39.  Ali  sensitivity 
analyses  tvere  performed  for  a period  of  100  years  using  the  AIT  sea  level  rise  scenario. 

The  hydraulic  conduclivi^  of  the  models  was  modified  for  a homogeneous  aquifer. 
Similarly,  tests  for  heterogeneity  were  performed  by  changing  the  oricnlalion  of  the 
individual  layers.  Layers  of  varying  hydraulic  conductivity  were  first  arranged 
horizontally  and  then  vertically  (see  Figure  7-kO). 


The  shift  in  the  imerfsce  was  plotied  for  seleci  values  of  horiionlal  hydraulic 
coftducliviry  Ki,.  The  results  rfiow  lha:  increase  in  Ks  has  moderate  effect  on  the 
simulnied  concentrations  and  heads.  Increasing  the  value  of  Ks  from  90  ni/day  to  9000 
m/day  resulted  in  shift  ofthe  0.1  isochlor  by  about  3500  m inland  (see  Figure  7-41). 
Similarly  reducing  the  values  of  K<  shifts  the  interface  towards  the  seaside  boundary. 


is  observed  in  Figure  7-41.  Forhi^er  vaiuesofisochiors.  especiaily  ciose  io0.95.  Uie 
change  seems  to  be  insignificam  as  the  isochlors  are  heid  at  the  coast.  Thus,  lower  vaitics 


of  Kt,  will  result  In  a i 


' for  the  same  set  of  aquifer  parameters  ortd 


effect  on  the  simuialed  concentration  values  (see  Figure  7>42).  The  isochlors  ^iAed  very 
liDle  for  a change  in  values  of  K.  ranging  from  I m/day  to  90  m/day. 


Sensitivity  to  aquifer  heterogeneity  was  also  tested  for  horizontal  and  vertical 
layers  of  varying  hydraulic  conductivities  as  described  in  Figure  7-40.  As  observed  in 
Figure  7-42,  the  vertically  layered  heterogeneous  aquifer  exhibited  appreciable  results  I 


comparably  well,  although  h was  not  as  significant  as  the  vertical  layering,  llie  practical 
significance  of  creating  physical  barriers  eilher  from  impervious  materials  or  layers  of 
very  low  hydraulic  cortductivlty  is  a method  being  investigated  by  researchers.  The 
creation  of  physical  barriers  helps  to  mitigale  the  impacts  of  seawater  intrusion  by 
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Figure  7-43.  Sensitivity  piots  for  aquifer  heterogeneity. 

Sensitivity  to  flux  boundaty  was  tested  by  modifying  the  specified  flux  boundaty 
on  the  [eft  hand  side  of  the  mode!.  As  observed  in  Figure  7-44,  increasing  the  magnitude 
orfiux  resulted  in  shifting  of  the  isochtors  seaward.  Fora  flux  of  IS  mVday/m  of 
coastline  the  0.1  isochlor  was  at  a distance  ofaboutdiOO  m while  increasing  the  flux 
magnitude  by  from  IS  to  30  mVday  resulted  in  shift  of  the  0.1  isochlor  by  about  i.OOO  m 
towards  the  const.  Decreasing  the  flux  magnitude  hod  the  opposite  effect  on  the  isochlors 
and  they  shifted  inland.  The  isochlors  less  than  the  0.6  are  observed  to  be  more  sensitive 
while  the  ones  close  to  0.9S  seem  to  be  locked  at  the  coast.  Compaietively  the  influence 
of  increasing  flux  is  a similar  to  reducing  horizontal  hydraulic  conductivity  and  vice 


significance  oFa  hydraulic  barrier  in  mitigating  the  impact  of  saltwater  intrusion.  Ifan 


possible  to  slowdown  the  process  of  aquifer  salinizalion.  In  it  is  possible  to  reverse 
the  process  of  salinization  by  creating  a hydraulic  blKk  (Tiruneb  and  Mots,  2003) 


Longitudinal  dispersivity  ql  was  found  to  influence  the  salinity  as  its  values  were 
increased.  The  sensitivity  exhibited  a significant  change  as  the  values  of  oc  approached 
the  rule  of  thumb  values  of  about  l/IOofihe  length  ofthecross-seciion.  Withan 


increase  in  ol.  the  salinity  values  tend  to  spread  out  farther  inland,  especially  within  the 
l/IO  length  range.  Transverse  dispersivity  ot  on  the  other  hand  does  not  seem  to  affect 
the  salinity  values  significantly.  With  increase  in  dispersivity.  the  isochlors  tend  to  be 
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more  straight  and  vertical.  Also,  when  the  dispersivity  value  was  kept  very  low,  the 


Two  types  of  sensitivity  analysis  were  performed  to  investigate  sensitivity  to 
recharge.  The  first  was  uniform  recharge  rate  applied  throughout  the  simulation  period, 
but  increased  for  different  tests.  The  second  was  variable  recharge  rate  that  ollemaied 
between  wet  and  dry  decades  during  the  entire  simuJation.  For  the  first  type  of  sensitivity 
analysis,  a sustained  uniform  rate  of  recharge  resulted  in  a shift  of  the  isochiors  seaward 
for  increasing  rates  of  recharge  (see  Figure  7-45).  For  the  second  type  of  analysis,  which 
Is  more  realistic  compared  to  the  first,  no  appreciable  change  in  the  position  of  isochiors 
was  noticed. 

It  is  understood  that  pumping  plays  a very  significant  role  in  the  dynamics  of 
saltwater  intrusion.  In  order  to  quantitatively  analyze  the  process  the  same  aquifer  was 
subject  to  varying  pumping  stresses.  Five  pumping  wells  were  used  to  pump  a total  of 
50. 500  and  5.000  m^/day  in  three  sensitivity  tests.  As  observed  in  Figure  7-46, 
increasing  the  pumping  mie  to  500  mVday  resulted  in  shift  of  the  0.1  isochlorby  less  than 
500  m,  while  pumping  at  5,000  mVday  resulted  in  complete  intrusion  of  the  aquifer. 
During  the  pumping  sensidvity  analysis,  the  flux  at  the  left  hand  boundary  was  kept  at  1 5 
mVday/m  of  coastline,  which  is  3,000  mVdoy  of  flow  for  the  200  m cell  width  used. 
Accelerated  intrusion  was  observed  when  the  amount  of  pumped  water  exceeded  the 
magnitude  of  flux.  This  type  of  sensitivi^  analysis  in  addition  to  other  oplimol 
management  methods  helps  to  determine  the  amount  of  water  that  could  safely  be 


withdrawn  for  use. 
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Figure  7*46.  Sensiiivity  plols  for  pumping. 


governing  flow  and  adveciion-dispersion  equations.  Experience  has  shown  that  soiulions 
obtained  (join  traditional  finite  difference  methods  quite  often  exhibit  numcricai  damping 
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or  spurious  oscillBUons.  This  nuiDerical  dispersion  couid  be  broughl  about  by  more  than 
one  factor.  The  widely  held  opinion  is  that  the  munerical  dispersion  is  caused  by  the 
inadequate  approximation  of  the  adveclion  terms  in  the  governing  transport  equation 
(Bear  and  Vemiijt.  1937;  Domenico  and  Schwartz,  1999;  Zheng  and  Bcnnei.  1995). 
Them  am  also  some  researchers  who  argue  that  the  problem  can  also  arise  due  to 
inaccurate  approximation  of  the  dispersive  fluxes,  more  specifically,  traditional  finite 
difrerence  methods  like  the  ones  used  in  MT3D.  SEAWAT,  RT3D  can  lead  to  significant 
unphysicai  oscillations  and  negative  concentrations  when  dispersion  is  suongly 
anisotropic  and  the  principal  direction  of  anisotropy  deviates  significantly  from  the  grid 
orientation  (Bear,  1972;  Bear  and  Bachniat,  1967,  Domenico  and  Schwartz  1999). 

Numerical  solutions  of  both  the  groundwater  flow  and  transport  equations  are 
bo.scd  on  a number  of  solution  techniques.  Modeling  techniques  used  to  solve  the  flow 
and  transport  problems  in  the  past  include  analog  models,  finite  difterencc.  finite 
elements,  analytic  elements,  method  of  characteristics  and  random  walk. 

The  finite  difference  method  implemented  in  bothMODFLOW  andMT3DMS  is 
used  in  SEAWAT  as  well.  The  finite  differenee  approximation  quite  often  exhibits 
numerical  dispersion,  which  is  caused  by  truncation  of  the  terms  that  approximate  the 
flow  in  the  formulation.  In  order  to  ensum  numerical  stability  and  accuracy  of  the 


adveclion-dispersion  solution  certain  stability  analysis  criteria  need  to  be  satisfied.  One 
such  criterion  Involves  performing  eigenvalue  analysis. 

In  order  to  successfiilly  complete  a simulation,  it  is  necessary  to  spend  lime 
experimenting  with  the  various  solution  techniques  and  alter  parameters  to  encourage 
convergence  and  minimiae  numerical  dispeision.  The  general  set  of  conditions  for 
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convergence  and  minimal  numerical  dispersion  are  sel  in  terms  of  flow  velocity,  flow 
lime  step,  and  aquifer  spatial  discretization.  TliediscreiiMUon  criteria  used  In  SEAWA'f 
are  based  on  that  of  MT3DMS.  The  stability  consirauils  and  accumey  requirements  for 
the  transport  of  conservative  species  with  the  explicit  Unite-difference  scheme  as  (Zheng 
and  Wang,  1908); 

adveclion,  , (7-3) 

dispersion,  (7-4) 

ix’  iy’  Ar' 

sink/source,  A/  S r^,  (7-5) 

k.| 

where  Jr  is  the  length  of  the  lime  step,  and  Ji,  Jy,  and  Je  represent  the  model  cell 
dimensions  in  the  x,y,  and  z directions.  The  stability  criteria  implemented  in  MT3DMS 
are  given  in  equations  7-3  to  7-5.  The  inherent  formulation  of  the  explicit  finite 
difference  scheme  allows  the  program  to  calculate  the  length  of  transport  lime  step.  The 
stability  constraint  for  adveclion  requires  flow  velocities  in  order  to  calculate  transport 
step  lengths. 

If  the  implicit  iterative  procedure  called  the  GCC  (generalized  conjugate  gradient) 
solver  is  used,  then  the  user  has  the  option  of  specifying  transport  time  step  lengths. 
Although  the  implicit  solver  may  reduce  the  number  of  transport  steps  required  fora 
simulation,  there  will  be  limitations  doc  to  convergence  and  accuracy  requiremenis. 
Generally  the  user  has  to  specify  smaller  lime  steps  in  onier  lo  encourage  convergence. 
Using  the  implicit  scheme  also  requires  observing  Ihc  stability  criteria  set  in  terms  of  the 
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dimensionless  grid  Peclet  number  end  Counmt  number  which  are  given  in  equations  7*6 
and  7-7 


C.  = 


PAT, 


(7-7) 


where  V is  flow  velocity  and  is  dispersion  coeiTicient.  [n  advection  dominated 
problems,  the  Peclet  mimberhas  to  be  kept  sufficiently  higher  to  ensure  stability, 

The  solution  schemes  provided  in  MT3DMS  include  Finite  DifFcrence  (FD). 
Method  of  Chamcteristics  (MOC),  Modified  Method  of  Characteristics  (MMOC),  Hybrid 
Method  ofChamcleristicsfHMOC),  and  Total  Variation  Diminishing  (TVD).  Generally 
the  FD  technique  is  faster  to  run  and  mote  vulnerable  to  convergence  problems  while  the 
MOC,  MMOC.  and  HMOC  have  slightly  better  sensitivity  to  numerical  convergence  and 
take  longer  time  to  rtin.  The  uldmate  method,  TVD  ensures  convergence  in  many  cases 
but  takes  a considerably  longer  lime  to  run  compared  toothers.  Also  the  FD  method 
exhibits  more  numerical  dispersion  compared  to  MOC  and  TVD.  although  the  TVD  is 
essentially  a higher  order  FD  technique  based  on  the  premise  of  a diminishing  sum  of 
concentration  differences  between  adjacent  nodes.  Sensitivi^  analysis  was  performed 


for  the  Coconut  Grove  model  using  all  the  solution  techniques.  The  results  are  presented 
in  figure  7.44  and  7-45. 
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Figure  7*4S,  Sensilivity  of  concenlrvion  profile  ro  solution  leciinique. 
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ApplicatiDD  of  the  Genetic  Algoritbms  Based  Opilinizatiun  Model 

The  oplimizaiion  modeling  was  tested  on  a surHcial  unconrined  aquifer  of  extent 
3000m  in  both  width  and  length.  The  aquifer  has  three  layers  and  a total  depth  of  60  in. 
The  model  was  discretized  into  30  columns  by  30  rows  of  each  100  m cell  size  (see 
Figure  T-dO).  The  seaside  boundary  is  set  to  have  a lime  variant  constant  head  and 
seawater  concenualioQ  of  35  kg/m^  The  intrinsic  aquifer  parameters  include  a horizontal 
hydraulic  conductivity  of  1000  m/day,  verticoi  hydraulic  conductivity  of  lOm/day, 
porosity  of  0J5,  longitudinal  dispeisivity  of  300  m,  transverse  dispersivity  of30mand 
vertical  disfieisivity  of  3 m.  A uniform  recharge  rate  of  1320  mm/year  was  applied. 

TwoEctsof  simulations  were  performed  by  operating  four  wells  at  a time.  The 
optimization  well  on/off  slalus  was  utilized  to  turn  four  wells  off  at  a given  time  and 
perform  the  optimization.  In  the  second  set  of  optimization  the  other  set  was  turned  on 
while  the  first  four  were  turned  off.  The  pumping  wells  were  spaced  1 000  m apart 
horizontally  in  order  to  minimize  the  chance  of  well  interference.  In  order  to  satisfy  the 
well  interference  criteria  wells  1 to  4 were  turned  on  initially.  During  the  second 
optimization  run.  wells  5 to  8 were  turned  on.  The  model  does  not  support  a moving  well 
option,  bul  this  approach  makes  it  possible  to  perform  opiimizadon  where  the  spatial 
location  of  the  well  plays  a role  in  determining  Ihc  optimal  strategy,  in  all  the 
optimization  runs,  pumping  wells  1. 2, 5.  and  6 pump  from  layer  2.  while  wells  3,4, 7, 
and  8 pump  from  layers  3 and  4. 


Figure  7>49.  Model  configuration  and  location  of  pumping  wella. 

The  management  period  considered  in  this  model  is  20  years  and  initial  simulations 
were  performed  to  see  the  level  of  salinity  in  each  well  and  generate  the  required  inputs 
for  the  optimization  model.  The  initial  pumping  and  salinity  data  are  presented  in  Tahie 


Well  ID 

Distonee  from  sea 

Pumping  rote 
tmVdavt 

Observed  salinity  range 

1800 

V 

2000 

2400 

1 

2000 

2800 

4 

1000 

3000 

Well  5 

1800 

Well  6 

2500 

Well? 

2500 

2800 

0 

* 

1500 

3000 

0-2.4 

concentration  values  had  to  be  set.  The  minimum  concentration  value  Is  set  to  be  0 while 


I value  is  set  to  be  O.S  kgfm^  The  maximum  value  is  based  on 
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the  requiremenU  for  polsbie  water  supply.  Also  the  minimum  pumping  rate  was  set  to  be 
0 while  the  maximum  pumping  rate  was  set  lo  be  ICT*  mVtlay.  In  order  to  estimate  the 
global  maximum  pumping  value,  the  net  recharge  was  calculated  and  divided  for  the 
existing  four  wells.  The  resulting  figure  was  higher  than  IxIO'mVd^  so  it  wasassumed 
that  a Ixl0‘  mVday  maximum  pumping  was  reasonable  to  use.  Perhaps  in  other 
situations  this  figure  could  be  set  depending  on  additional  factors  such  as  forccasred 
demand. 

Optimal  Pumping  Strategy 

The  optimisations  performed  for  the  two  layouts  of  wells  required  an  average  of 
700  generations  with  a crossover  probability  of  0.85  and  the  option  ofno  elitian.  The 
results  of  optimal  pumpbg  strategy  are  presented  in  Table  7-10. 


Table  7-10.  Results  of  optimal  pumping  strategy  aiui  observed  aquifer  salinity. 


Well  ID 

Initial  pumping 
(mVday) 

Optimal  pumping  rates 
(m^/day) 

Range  of  Salinl^ 
(kg/ra‘) 

Layout  1 

Layout  2 

Layoirl  1 

Layout  2 

Well  1 

1800 

800 

0-0.5 

Well  2 

2400 

3000 

0 

Well} 

2800 

3000 

0 

Well  4 

3000 

200 

0-057 

Well  5 

1800 

1000 

0-03 

Well  6 

2400 

2400 

0 

Well  7 

2800 

2800 

0 

Well! 

3000 

1000 

0.55 

The  distribution  of  salinity  isochlots  and  hydraulic  heads  for  pre  and  post 
oplimiaalion  periods  are  presented  in  Figures  7-50  and  7-51  for  layout  1.  end  Figums  7- 
52  and  7-53  for  layout  2. 


Figure  7‘50.  Dislribulion  of  salinity  isochlors  for  the  current  (a)  and  optimal  (b)  pumping 
strategy  at  the  end  of  20  years  management  period  for  Layout  I. 


Figure  7-5 1 . Distribution  of  groundwater  head  for  the  current  (a)  and  optimal  <b) 

pumping  strategy  at  the  end  of  20  years  management  p^iod  for  Layout  I. 

From  the  concentiarion  plots  presented  in  Figure  7-50  a and  b,  it  is  possible  to  see 
the  inipuclof  optimized  withdrawal  on  the  process  of  saltwater  intrusion.  If  layout  I is  to 


be  adopted,  the  current  pumping  rale  of  10,000  ra'/day  has  to  be  reduced  to  7.000  mVday, 


'This  results  in  a slower  rate  of  saltwater  intrusion  by  mainlainingtheO.S  ppm  isochlorat 
a mtutimum  distance  of  about  1,500  m Inland.  Layout  2 on  the  other  hand  gave  an 


optimal  pumping  rate  of  7,200  m^/day  for  the  * 
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moniloring/pLimping  points.  The  plot  of  changes  in  concentration  for  both  layouts  is 
presented  in  Figure  7*54. 


From  the  concenlraiicm  plots  presented  in  Figure  7'50a  and  b,  it  is  possible  to  see 
the  impactoroptimized  withdrawsi  on  the  process  of  saltwater  intrusion.  IflayoutI  is  to 
be  adopted,  the  current  pumping  rale  of  10,000  m^fday  has  to  be  reduced  to  7,000  mayday. 
This  results  in  a slower  rate  of  saltwater  intrusion  by  maintaining  theO.S  ppm  isochlorat 
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a maximum  distance  of  about  1 ,500  m inland.  Layout  2 on  the  other  hand  gave  an 
optimal  pumping  rate  of  7,200  mVday  for  the  same  criteria  of  lower  salinity  at  the 
motiitoringfpumptng  points.  The  plot  of  changes  in  concentration  for  both  layouts  is 
presented  in  Figure  7-54. 


Figure  7-54.  Plot  of  change  in  concentration  and  pumping  for  Layoutl  and  Layout  2. 

The  plot  of  changes  in  pumping  and  monitored  concentration  are  plotted  in  Figure 
7-54,  exhibiting  an  utcreasing  trend.  It  is  also  possible  to  conclude  that  location  of  the 
wells  plays  a significant  role  in  the  solution  of  the  optimization  problem.  The  pumping 
values  prcscoled  ate  by  no  means  global  optimum  values,  but  rather  close  to  what  could 
be  global  maximum  given  the  nature  of  the  algorithm  used  in  the  solution. 

The  salt  mass  extracted  through  pumping  was  one  of  the  constraining  variables.  As 
presented  In  Table  7- 1 1 . the  sail  mass  extracted  during  post  optimization  simulation  Is 
considerably  lower  thim  the  pre.optimiznlion  simulation.  A tradeolT curve  was  also 


plotted  for  both  layouts  to  see  the  variation  between  salt  mass  and  pumped  volume  for 
both  (see  Figure  7-S5). 


Salt  mass  extracted  and  volume  o 


Silt  Miss  EalractedthB) 


Volume  Pumped  lm^> 


As  observed  from  the  slopes  of  the  linear  equations  in  Figure  7*5S.  layout  2 
provides  a smaller  amount  of  salt  for  the  volume  of  water  pumped.  The  equations 
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possible  a well  informed  decision. 


130 


h is  also  important  to  incorporate  non-quantifiabie  socio'ocononiic  factors  in  the 
developmem  of  the  decision  support  system.  This  could  be  achieved  by  involving  the 
various  stakeholders  in  the  management  of  the  water  resources  of  the  area. 


CHAPTERS 

SUMMARY  AND  CONCLUSIONS 
Ad  EvalualioD  of  (Uc  Roscarcb  PresoDted  Ibis  DisserialioD 
This  research  is  the  result  of  about  three  and  half  years  ofsvork  that  had  the  purpose 
of  esl^lishiRg  a sound  understanding  of  the  issues  of  climate  change  induced  sea  level 
rise,  saltwater  intrusion  in  coaslal  aquifers,  and  integrated  management  strategy  for 
coastal  aquifers  on  the  basis  of  optimal  withdrawal  of  freshwater. 

On  Ihe  basis  ofihe  objectives,  a GIS  based  tool  that  can  be  used  to  investigate  the 
impacts  of  climate  change  induced  sea  level  rise  was  developed.  This  approach  also 
leaves  room  for  further  expansion  of  research  to  include  integration  of  other  components 
ofeoastaJ  hydrologic  processes.  The  three-dimensional,  variable  density,  groundwater 
now  modeling  code  SEAWAT  was  modified  to  incorporate  a sea  level  rise  component 
that  dynamically  updates  the  coaslal  boundary  to  simulate  the  rise  in  sea  level  as  part  of 
Ihe  hydrostatic  time  variant  specified  head  sea  boundary.  An  oplimizaiian  code  also  was 
developed  on  Ihe  basis  ofihe  mechanics  of genetic  algorilhms.  The  oplimi^aIion  code 
uses  the  outputs  of  the  simulation  model  to  compute  optimal  pumping  values. 

The  Ihree  simulalion-oplimization  modeling  tools  me  combined  in  a shell  program 
that  enables  seamless  running  of  each  component  within  a uscr-lriendly  graphical  user 
Interface.  In  addition,  a bamework  that  could  be  used  to  develop  a decision  support 
system  for  optimal  management  of  coastal  aquifers  subject  to  climate  change  impacts 
also  was  developed.  Detailed  summmy  and  conclusions  drawn  from  the  research  are 
presented  under  the  following  subheadings. 
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Global  ClimaK  Cbangc^ea  Levtl  Rise.and  Soawaler  Inlrusion 
There  is  a compelling  evidence  suggesting  the  signilicance  of  anlhropogcnie  foelors 
in  global  climale  and  change.  One  of  the  primary  manifestoiiona  of  anthropogenic 
climaicchangeis  the  rise  in  the  global  sea  level.  The  sea  level  Is  expecled  lo  rise  mainly 
due  10  global  warning  and  Ihc  subsequent  melting  of  ice  caps  and  thermal  expansion. 

The  most  recent  climate  change  impacts  assessment  (IPCC,  2001)  concluded  that  the 
global  sea  level  is  projected  to  rise  by  0.09  to  0.88  meters  between  1990  and  2100  for  ihe 
lull  range  ofSRES  scenarios.  The  projcclions  made  by  the  IPCC  are  slightly  lower  lhan 
the  earlier  forecasts  made  on  the  basis  of  the  IS92  marker  scenarios,  which  projected  the 
sea  level  lo  rise  globally  by  0. 1 3 lo  0.94  meters. 

Climatic  anomalies  such  us  El  Nido  arc  expected  to  continue  with  little  or  no 
change  in  their  amplitude.  The  impacts  of  El  Niito  arc  primarily  extremes  of  dry  weather 
and  heavy  rainfall  depending  on  the  geographic  location  of  the  area.  Extremely  dry 
weather  will  have  significant  impacts  on  Ihe  estimation  of  actual  and  potential 
evapotranspiralion  while  the  high  rainfall  will  pose  a serious  flooding  threat.  Global  sea 
level  rise,  on  the  other  hand,  has  a delayed  response  time  to  such  an  intense  spell  of  wet 
period,  and  may  not  be  signiEcanlly  affected  in  the  short  run.  However,  the  development 
of  a management  model  for  coastal  aquifers  willi  a management  period  of  decadal  to 
century  scale  has  to  consider  the  impacts  of  El  Niho  and  other  climatic  anomalies  in  the 
form  of  changes  in  precipitation  and  evaporation. 

Sea  level  rise  is  likely  lo  cause  accelerated  seawater  intrusion  on  a scale  of  century 
or  more.  When  the  aquifer  is  subject  to  additional  stiesses  due  to  pumping  and  reduced 
recharge,  the  rate  of  intrusion  Is  accelerated  even  foster.  The  influence  of  sea  level  rise 
and  freewater  withdrawal  on  the  process  of  seawater  intrusion  Is  dependent  on  a number 
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orfaclors.  The  faclors  include  aquifer  hydraulic  properties,  aquifer  geometry,  the 
location  of  pumping  wells,  and  the  magnitude  of  pumping. 

Based  on  the  analysis  conducted  for  the  surficial  aquifer  system  of  south  Florida, 
climate  change  and  sea  level  rise  pose  a serious  threat  to  the  ecosystem.  Within  the  neat 
100  years  or  more,  the  primary  forms  of  changes  in  the  ecosystem  are  likely  to  be 
accelerated  translation  of  the  shoreline.  Vulnerability  assessmeou  for  the  three  counties 
considered  in  this  study  indicate  that  at  a scale  of  a century  the  relative  sea  level  rise 
esiimuiB  for  the  south  Florida  coast  will  be  about  O.S  m.  For  the  entire  range  of  IPCC 
scenarios  considered  for  select  tidal  slations.  the  relotive  sea  level  rise  was  estimated  to 
be  between  O.d  to  0.6  m.  Depending  on  Ihe  slope  of  the  local  area  considered,  the 
corresponding  shoreline  translation  for  Palm  Beach  County  could  range  from  a few 
meters  to  about  200  meters  or  more.  Similarly,  for  Broward  County,  the  estimated 
shoreline  translation  could  be  as  high  as  1,000  m.  For  coastal  areasofMiami-Dode 
County,  the  shoreline  Uanslalion  could  be  more  lhan  100  m. 

This  study  serves  as  indicator  of  the  vulnerability  of  coastal  areas  to  sea  level  rise. 
Detailed  assessments  should  include  mapping  of  existing  infrastructure,  settlement 
putlems,  natural  habitat,  wildlife,  and  water  resources.  The  study  should  also  be 
performed  al  a smaller  scale  lo  enhance  the  visibility  of  the  different  GIS  dala  layers  to  be 
used  in  Ihe  analysis.  A delailed  study  should  also  consider  that  a range  of  scenarios  and 
selection  has  to  be  made  on  the  basis  of  likelihood  of  occurrence. 

The  analysis  of  available  topographic  data  for  Goksu  Delta  in  Turkey  indicated  that 
sea  level  rise  is  a factor  that  needs  lo  be  considered  in  future  studies.  Detailed  sea  level 
rise  and  shoreline  Uanslalion  analyses  were  not  conducted  for  the  Goksu  Delta  due  to 


unavailability  of  a DEM  for  the  area.  As  a result  it  is  not  possible  to  make  specific 
conclusions. 

Sea  level  rise  is  also  noted  to  haves  significant  impact  on  the  movement  of  the 
transition  zone.  For  the  cross-section  model  studied  in  south  Florida,  portion  of  the 
aquifer  that  catends  from  3.000  to  4.200  meters  from  the  coast  is  expected  to  undergo 
significant  changes  in  ambient  aquifer  salinity  due  to  sea  level  rise.  In  general  there  will 
be  increased  levels  of  salinity  and  shili  of  the  transition  zone.  The  general  trend  of  flow 
is  not  expecled  to  show  significant  changes.  The  groundwater  flow  regime  will  be 
afTecied  in  terms  of  slightly  higher  levels  of  hydraulic  head  and  a more  dynamic 
exchange  of  ffeshwater-saltwater  near  the  coast. 

The  surlicial  aquifer  system  of  south  Florida  is  not  likely  to  attain  a stale  of 
dynamic  equilibrium  withm  the  next  millennium.  Pumping  streasos  and  the  system  of 
canals  that  exist  in  the  area  are  the  primary  factors  that  make  it  difficult  to  attain  dynamic 
equilibrium.  Detailed  analyses  that  consider  pumping  withdrawals  and  ihe  full  range  of 
sea  level  rise  scenarios  are  expected  to  give  a more  practical  management  scenario  for  the 
next  century. 

The  study  conducted  on  the  Goksu  Delta  aquifer  indicated  that  the  change  in 
aquifer  salinity  is  expected  to  cover  a cross-section  Ihai  is  about  1.000  meters  wide  on 
average  in  the  X-X  section  (north-south  orientation).  The  Y-Y  section  feast-west 
oricnialion)  exhibited  a much  wider  zone  of  influence  due  to  the  existence  of  a seaside  on 
both  ends.  A ship  ofaquifer  that  has  a width  of  about  10,000  meters  and  a depth  that 
extends  from  200  to  400  meters  is  expected  to  undergo  significant  changes  In  salinity  due 
to  sea  level  rise.  The  significance  of  geometry  is  also  very  noticeable  in  the  Goksu  Della 


i observed  by  Ihe  shape  of  the  transilion  of  the  zone.  Conclusions  similar  to  the  one  Tor 


thesurriciol  aquirer  eioss-scction  orsouth  Florida  could  be  made  Tor  Ihe  Goksu  Delta 
aquifer. 


Generic  algorithm  based  search  and  optimization  methods  provide  a very  robust 
alternative  to  thccaisting  optimization  methods.  The  nonlmeur  nature  of  most  water 
resources  problems  makes  non  gradient-based  search  and  optimization  methods  more 
appealing.  The  selection  of  parameters  and  seeding  of  Ihe  random  number  generation 
scheme  require  careful  understanding  and  experimentation.  The  stopping  criterion  needs 
to  be  dchned  on  Ihe  basis  of  a more  cfficiem  method.  Generally,  genetic  algorithm  bused 
optimization  requires  longer  piocossing  time.  Therefore,  access  to  faster  processors 
makes  the  process  run  smoothly. 

The  research  presented  in  this  dissertation  is  based  on  exlemaJ  coupling.  However, 
the  process  of  simulation-optimization  runs  much  better  if  the  two  processes  are  coupled 
intemally.  With  internal  coupling,  the  simulation  results  arc  transferred  dynamically  to 
Ihe  optimization  program  and  vice-versa. 

Modeling  IssuesKkb^ives  and  Seale  Coosideratioos 
There  are  a number  of  approaches  to  modeling  seawater  intrusion.  One- 
dimensional.  steady  slate,  analytical  expressions  such  as  tlie  Ghyben-Herzberg 
approximation,  are  simple  to  use  and  do  not  require  much  data.  They  also  have  link 
computational  requirements.  The  premise  of  this  analysis  is  the  existence  of  steady  state; 
however,  that  is  not  quite  the  case  os  seen  in  the  sensitivity  analyses.  These  types  of 
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analyses  involve  limiting  assumptions  and  could  be  performed  in  a relatively  short  penod 

Two-dimensional,  transient  mimcHcal  models  on  the  other  hand  can  handle  more 
complex  situations  and  provide  a more  realislic  representation  of  the  processes  and  can 
account  for  variability  in  physical  situations.  They  have  the  ability  to  forecast  and 
develop  scenarios  by  considering  temporal  variability.  Two-dimensional  models  require 
computer  modeling  and  involve  moderate  to  targe  computational  efToit.  There  are  some 
analytical  approximalions  for  two-dimensional  problems.  Being  numerical 
representations,  they  could  have  numerical  stability  problems  and  require  proper 
undcrslanding  of  solution  methods  and  techniques.  The  processing  time  is  relatively 
grealer  lhan  for  the  one-dimensional  models,  generally  ranging  from  a few  hours  to 
several  days  depending  on  the  processor  speed. 

Three-dimensional  numerical  models  provide  a much  more  accurale  representation 
of  the  aquifer  system  to  be  modeled-  Three-dimensional  models  require  intensive 
computational  efforts  and  considerably  longer  simulation  lime.  The  results  arc  much 
more  sensitive  to  temporal  and  spatial  discretization  and  solution  techniques  and  more 
vulnerable  to  numerical  dispersion.  As  ills  normally  die  case  with  variable  density  three- 
dimensional  models  of  moderate  extent,  the  computing  needs  are  also  coupled  with  input 
data  and  output  file  storage  needs.  They  generally  require  a thorough  understanding  of 
the  numerical  methods  involved  in  order  to  obtain  a mimerically  stable  solution.  The 
results  are  more  useful  and  accurale  in  predicting  the  future  and  developing  scenahos. 

Sharp  interface  models  are  useful  alternatives  that  provide  insight  into  the  flow 
regime  and  characteristics  in  coastal  aqaifers.  However,  the  development  of  sharp 
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imcrfoce  models  needs  lo  be  perfcmned  with  care  and  understanding  of  the  characteristics 
of  both  the  area  that  needs  to  be  modeled  and  the  underlying  assumptions  of  the  sharp 
interface  approach.  Perhaps  the  two  most  significant  factors  that  need  to  be  considered 
are  the  thickness  of  the  interface  and  equilibrium  in  the  aquifer.  Generally,  sharp 
interface  modeling  should  not  be  applied  if  the  interface  is  thick  compared  to  the  depth  of 
the  aquifer.  The  second  signiheant  factor  is  the  assumption  ofa  hydrostatic  equilibrium, 
which  in  effect  is  not  easy  to  achieve  as  aquifers  are  generally  m a dynamic  equilibrium 

Experience  has  shown  that  most  coastal  aquifers  have  a significantly  broad 
uansition  zone.  Thus,  using  the  sharp  interface  approach  does  not  provide  an  accurate 
representation  of  the  flow  process.  In  silualions  where  there  is  a freshwater  lens  or  a 
transition  zone  that  extends  to  a few  meters,  then  the  sharp  interhice  could  be  used  by 
assuming  the  Ghyben-Herzberg  relaiionship. 

Management  models  forcoaslal  aquifers  can  have  varying  lime  scales.  If  the 
management  lime  frame  is  not  more  than  a few  decades,  then  the  impacts  of  climate 
change  can  be  approximated  by  the  variability  in  reohaige  and  evapotranspiratioin.  If  on 
the  other  hand,  the  lime  frame  is  on  a scale  of  a century  then  the  impacts  of  sea  level  rise 
need  to  be  incoiporaied  in  addition  to  variability  in  recharge  and  evapotnmspi  ration. 

Water  resources  management  has  to  follow  an  integrated  approach  that  considers 
both  surface  and  groundwater  sources  as  an  integral  unit  Depending  on  the  scale  of 
management  time,  the  cumulative  impact  of  climate  change  and  consumptive  use  has  to 
be  incorporated  in  the  development  of  a regional  water  resources  analysis  model. 
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Development  of  a framework  for  management  of  coastal  water  resources  has  to  be 
focused  on  optimal  withdrawal  of  h'eshwaler  from  both  surface  artd  groundwater.  Since 
the  problem  of  seawater  intrusion  is  evident  in  most  coastal  aquifers,  the  manogement 
practice  needs  to  have  a component  of  impact  mitigation.  This  could  be  accomplished  by 
utilising  a number  oflechniques  such  as  a hydraulic  barrier,  physical  barrier,  and  surface 
recharge.  Development  of  a regional  model  helps  to  formulate  a decision  support  tool 
that  enables  management  of  water  resources  based  on  a sound  policy  that  ensures 
sustainable  use  of  water  resources.  A flowchart  for  a decision  support  system  that  could 
be  usml  to  formulate  a water  resources  management  policy  for  coa-stal  areas  is  presented 
in  Figure  8-1. 
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APPENDIX  A 

DEVELOPMENT  OF  INTEGRATED  MODELING  ENVIRONMENT 
CSCAMM:  CUmaie  Sensitive  Consul  Aquifer  MaDagcment  Model 

CSCAMM  Is  a shell  prognun  wrilien  in  Visual  Basic  lo  serve  as  a Graphical  User 
Interface  (GUI)  and  faclliiaic  the  calling  ofthe  various  components  oflbe  model. 
CSCAMM  provides  the  user  the  ability  to  select  and  run  a specific  model  of  choice  by 
calling  the  execulables  specified  In  a selected  working  directory. 

CCGEN  Is  a CIS  tool  written  in  the  programming  language  Avenue  and  designed 
10  work  in  Arc  View,  therefore  the  user  has  to  make  sure  the  computer  used  for  the 
modeling  exercise  has  AreView  and  the  Spatial  Analyst  extension  is  activated.  The 
current  release  is  a fiilly  functional  program  lo  be  run  as  a project  In  Arc  View.  In  the 
future  CCGEN  will  be  recompiled  lo  be  an  AreView  extension,  thus  ensuring  case  in 
portability. 

SEAWATSLR  is  a Fonran  code  designed  to  enable  consideration  of  climate 
change  in  terms  of  sea  level  rise.  The  code  has  an  additional  rauiine  to  import  and 
analyze  sea  level  rise  data  and  modify  the  algorithm  for  Time  Varying  Constant  Head 
(CHD)  which  needs  to  be  defined  along  the  seaside  boundary.  During  the  modeling 
process  variability  in  recharge  and  evapotranspiralion  will  have  to  be  included  as  required 
for  the  area  of  study. 

GENOPT  is  a muHi-objeclivc.  multi-management  period  groundwater  flow 
opiimizaiicm  code  writlen  in  Foriraa  The  inputs  for  GENOPT  are  obtained  from  either 
SEA  WAT  or  SEAWATSLR. 
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INTERFACE  is  a Fortran  code  written  to  read  the  binary  oulpuls  ofSEAWAT  or 
SEAWATSLR  and  calculate  enviroiunental  groundwater  heads.  It  also  isolates  specified 
ieochlors  of  choice  and  enable  exporting  them  into  an  ASCII  file  which  could  easily  be 
imported  into  any  plotting  software  for  display. 

Software  and  Data  Requirements 

The  software  requirements  for  CSCAMM  are  primarily  of  two  types.  The  first  type 
of  software  is  a pre  and  post  processor  suited  for  MODFLOW  and  MT3DMS  based 
models.  In  this  research  the  pre  and  post  processing  was  accomplished  by  Groundwater 
Vistas,  a commercial  software  developed  by  Environmental  Simulations  Inc. 

The  second  type  of  software  is  the  Environmcniiil  Systems  Research  Institute 
(ESRI)  product  ArcVlew.  The  primary  advantage  of  AreView  is  its  availability  and  the 
ease  with  which  it  could  be  modified. 

Thedala  required  to  run  Ihe  various  components  of  CSCAMM  are  the  following: 

1 . Global  climate  change  scenario  to  be  selected  from  Ihe  IPCC  2001  documentation. 
This  ia  included  in  the  database  of  CSCAMM. 

2.  Digital  Elevation  data  to  be  obioined  from  appropriate  DEM  data  repository. 

3.  Local  sea  level  rise  trend  and  possibly  submergence  data  which  has  to  be  obtained 
from  NOAA  and/or  other  agencies. 

4.  Hydrngcologic  data  required  to  conslrucl  a variable  density  groundwater  flow 

5.  Water  demand  and  population  forecast  data  required  to  analyze  the  demand  during 
simulation  and  optimization. 

The  methodology  and  data  requiremems  are  outlined  in  figure  AA-l. 
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Process  Flowchart 


Computing  Requiracnenls 

The  computer  codes  used  in  this  research  are  all  designed  to  utilize  the  computing 
capacity  of  desktop  computers  in  a more  cmcienl  manner.  However,  due  to  the 


he  met  for  the  models  to  run  smoothly. 


processing  speed  for  the  simulations  to  run  eiricienUy.  In  addition  CCGEN  requires  a 
memory  of  ai  least  128  MB  in  order  to  be  able  to  run  a series  of  simulations.  Depending 
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or  Ihe  size  of  lie  raster  file  that  is  being  processed,  ihe  storage  capacity  of  Uie  disk  might 
sometimes  be  not  sufficient  resulting  in  a very  common  error  message  of  not  being  able 
to  create  grid  file  from  the  theme.  Unless  the  data  has  some  inherent  problems  this 
usually  comes  from  Insufficient  disk  space  or  not  the  theme  active  in  the  active 
document.  When  this  happens  ihe  ASCII  file  that  are  necessary  for  the  new  grid  files  to 
be  created  will  not  be  successfully  created  Ibus  giving  Ihe  error  message. 

The  user  has  Ihe  option  of  defining  the  analysis  side  depending  on  the  orientation 
of  Ihe  coast.  Since  the  translation  of  coast  is  calculated  from  the  sea  level  to  the  inland, 
the  user  has  to  specif  which  direction  during  the  eaecurion  of  the  program.  The  code  Is 
designed  to  be  Interactive  lo  allow  the  user  to  make  the  necessary  selections.  The 
simulation  using  CCQEN  also  creates  a file  named  “risejb*"  by  default.  This  file 
contains  the  necessary  information  required  lo  run  SEAWATSLR.  Upon  successful 
eseculion  ofthe  program  Ihe  intermediate  ASCII  files  will  bedeleled  thus  freeing  up 
space  that  could  unnecessarily  be  held  by  slope  and  elevation  ASCII  files. 


Figure  A-2.  Screen  capture  of  CCGEN  simulation. 
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Figure  A*  1 8.  Miami-Dade  coast  sea  level  Torecasl  for  the  year  2050. 
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